Methods and systems for updating a predicted location of an object in a multi-dimensional space

ABSTRACT

Embodiments of the present invention characterizing the uncertainty of the orbital state of an Earth-orbiting space object hereof using a Gauss von Mises probability density function defined on the n+1 dimensional cylindrical manifold    n × . Additionally, embodiments of the present invention can include transforming a Gauss von Mises distribution under a diffeomorphism and approximating the output as a Gauss von Mises distribution. Embodiments of the present invention can also include fusing a prior state represented by a Gauss von Mises distribution with an update report, wherein the update can be either another Gauss von Mises distribution of the same dimension as the prior or an observation related to the prior by a stochastic measurement model. A Gauss von Mises distribution can be calculated from a plurality of reports, wherein the reports are either Gauss von Mises distributions or observations related to the state space by a stochastic measurement model.

STATEMENT AS TO RIGHTS TO INVENTIONS MADE UNDER FEDERALLY SPONSOREDRESEARCH AND DEVELOPMENT

This invention was made with government support under Contract No.FA9550-12-C-0034 awarded by the Air Force Office of Scientific Researchto Numerica Corporation. The government has certain rights in theinvention.

FIELD OF THE INVENTION

The present invention relates generally to non-linear estimation andcomputerized techniques to characterize uncertainty on a cylindricalmanifold and, more particularly, to characterizing the state uncertaintyof Earth-orbiting space objects to support many functions in spacesituational awareness.

BACKGROUND OF THE INVENTION

Space Situational Awareness (SSA) includes knowledge of the near-Earthspace environment which can be accomplished through the tracking andidentification of Earth-orbiting space objects to protect space assetsand maintain awareness of potentially adversarial space deployments.Although current operational systems have performed well, future needswill far exceed current capabilities. With the instantiation of moreaccurate sensors and the increased probability of future collisionsbetween space objects, the potential number of newly discovered objectsis likely to increase by an order of magnitude within the next decade,thereby placing an ever increasing burden on current operationalsystems. Moving forward, the implementation of new, innovative,rigorous, robust, and autonomous methods for space object identificationand discrimination are required to enable the development andmaintenance of the present and future space catalog and to support theoverall SSA mission.

Fundamental to the success of the SSA mission is the rigorous inclusionof uncertainty in the space surveillance network (SSN). The propercharacterization of uncertainty, sometimes called covariance realism, isa common requirement to many SSA functions including tracking and dataassociation, resolution of uncorrelated tracks (UCTs), conjunctionanalysis and probability of collision, sensor resource management, andanomaly detection. While tracking environments, such as air and missiledefense, make extensive use of Gaussian and local linearity assumptionswithin algorithms for uncertainty management, space surveillance isinherently different due to long time gaps between updates, highmis-detection rates, non-linear and non-conservative dynamics, andnon-Gaussian phenomena. What is state-of-the-art and robust for air andmissile defense need not be applicable to SSA.

The field of sequential state estimation has much of its origins in thepioneering work of R. E. Kalman. Considered to be one of the most simpledynamic Bayesian networks, the Kalman filter updates a system staterecursively over time using incoming measurements and mathematicalprocess models. The basic Kalman filter assumes linearity in theunderlying dynamical and measurement models and that all error terms areGaussian. At the other end of the spectrum is the general Bayesiannon-linear filter which updates the full Probability Density Function(PDF) of the system recursively over time while permitting non-Gaussianerror terms and non-linear process models. Between the Kalman filter andthe general Bayesian framework are a host of sub-optimal methods forsequential filtering which have been developed over the past fifty yearsand tailored to specific applications. The most common extensions andgeneralizations of the Kalman filter include the extended Kalman filter(EKF) and the unscented Kalman filter (UKF) which both work onnon-linear systems. A feature of the latter is its “derivative-free”nature. Within the filter prediction step, the propagated state estimateand covariance are reconstructed from a deterministically chosen set of“sigma points” propagated through the full non-linear dynamics. Theequivalence between this reconstruction or “unscented transform” withGauss-Hermite quadrature has been established. Variations and extensionsof the UKF, including a more numerically stable “square-root” versionand the higher-order Gauss-Hermite filters, have been formulated. Inspace surveillance, the state or orbital uncertainty can be highlynon-Gaussian and filtering techniques beyond the EKF and UKF aresometimes required. Examples include Gaussian sum (mixture) filters,filters based on non-linear propagation of uncertainty using Taylorseries expansions of the solution flow, particle filters, and theprobability hypothesis density filter and its generalization, thecardinalized probability hypothesis density filter.

A drawback of many existing methods for non-linear filtering anduncertainty management, including the ones reviewed above, is theconstraint that the state space be defined on an n-dimensional Cartesianspace

^(n). Any statistically rigorous treatment of uncertainty uses PDFsdefined on the underlying manifold on which the system state is defined.In the space surveillance tracking problem, the system state is oftendefined with respect to orbital element coordinates. In thesecoordinates, five of the six elements are approximated as unboundedCartesian coordinates on

⁵ while the sixth element is an angular coordinate defined on the circle

with the angles θ and θ+2πk (where k is any integer) identified asequivalent (i.e., they describe the same location on the orbit). Thus,more rigorously, an orbital element state space is defined on thesix-dimensional cylinder

⁵×

. Indeed, the mistreatment of an angular coordinate as an unboundedCartesian coordinate can lead to many unexpected software faults andother dire consequences.

BRIEF SUMMARY OF THE INVENTION

Embodiments of the present invention overcome the disadvantages andlimitations of the prior art by providing methods and systems forproperly characterizing the uncertainty of a state defined on a n+1dimensional cylindrical manifold

^(n)×

, e.g., characterizing the uncertainty of the orbital state of anEarth-orbiting space object.

In order to provide a more statistically rigorous treatment ofuncertainty in the space surveillance tracking environment and to bettersupport the aforementioned SSA functions, embodiments of the presentinvention provide a new class of multivariate probability densityfunctions (PDFs), called Gauss von Mises (GVM) distributions, formulatedto more accurately characterize the uncertainty of a space object'sstate or orbit. Using the GVM distribution as input, extensions andimprovements can be made to tracking algorithms including the Bayesiannon-linear filter used for uncertainty propagation and data fusion,batch processing and orbit determination, multi-dimensional quadrature,and likelihood ratios and other scoring metrics used in dataassociation. An exemplary application of embodiments described herein isthe use of the GVM distribution within a sequential non-linear filterresulting in improved algorithms for uncertainty propagation and datafusion to support higher level SSA functions. Embodiments of the presentinvention provide a tractable implementation of the Bayesian non-linearfilter using the GVM distribution to model state or orbital uncertaintywhile providing significant improvements, both in terms of accuracy andcomputational expense, over the aforementioned methods.

In accordance with an embodiment of the present invention characterizingthe uncertainty of the orbital state of an Earth-orbiting space objecthereof represents the uncertainty using the Gauss von Mises probabilitydensity function defined on the n+1 dimensional cylindrical manifold

^(n)×

. Additionally, embodiments of the present invention can includetransforming a Gauss von Mises distribution under a diffeomorphism andapproximating the output as a Gauss von Mises distribution. Embodimentsof the present invention can also include fusing a prior staterepresented by a Gauss von Mises distribution with an update report,wherein the update can be either another Gauss von Mises distribution ofthe same dimension as the prior or an observation related to the priorby a stochastic measurement model. A Gauss von Mises distribution can becalculated from a plurality of reports, wherein the reports are eitherGauss von Mises distributions or observations related to the state spaceby a stochastic measurement model. Embodiments can include, but are notlimited to, providing a suite of decision-making tools, based on amethod for properly characterizing the uncertainty of the orbital stateof an Earth-orbiting space object as described herein, that allow anoperator to take an appropriate course of action such as maneuvering asatellite in order to avoid a potential collision or updating a catalogof space objects in order to improve the accuracy of the objects in thecatalog.

BRIEF DESCRIPTION OF THE DRAWINGS

The accompanying drawings, which are incorporated in and form a part ofthe specification, illustrate embodiments of the present invention and,together with the description, serve to explain the principles of theinvention. In the drawings:

FIG. 1 is a flowchart illustrating an exemplary process for transforminga Gauss von Mises distribution under a diffeomorphism and approximatingthe output as a Gauss von Mises distribution according to one embodimentof the present invention;

FIG. 2 is a flowchart illustrating an exemplary process for fusing aprior state represented by a Gauss von Mises distribution with an updatereport, wherein said update is another Gauss von Mises distribution ofthe same dimension as the prior according to one embodiment of thepresent invention;

FIG. 3 is a flowchart illustrating an exemplary process for fusing aprior state represented by a Gauss von Mises distribution with an updatereport, wherein said update is an observation related to the prior by astochastic measurement model according to one embodiment of the presentinvention;

FIG. 4 is a flowchart illustrating an exemplary process for generating aGauss von Mises distribution from a plurality of reports, wherein saidreports are Gauss von Mises distributions according to one embodiment ofthe present invention;

FIG. 5 is a flowchart illustrating an exemplary process for generating aGauss von Mises distribution from a plurality of reports, wherein saidreports are observations related to the state space by a stochasticmeasurement model according to one embodiment of the present invention;

FIGS. 6 a and 6 b are graphs illustrating the fusion of a priorprobability density function p₁(θ) with an update probability densityfunction p₂(θ) wherein (a) the angle θ is mis-represented as a Cartesianvariable and (b) the angle θ is represented correctly as a circularvariable according to one embodiment of the present invention;

FIG. 7 is a graph showing plots of the von Mises probability densityfunction with location parameter α=0 and different values of theconcentration parameter κ according to one embodiment of the presentinvention;

FIGS. 8 a and 8 b are graphs illustrating the setup for testing if apoint is a statistically significant realization of a Gauss von Misesdistribution with the shaded regions depicting those regions in whichthe null hypothesis is rejected according to one embodiment of thepresent invention;

FIGS. 9 a-9 f are graphs illustrating the uncertainty of a spaceobject's orbital state at different epochs t−t₀ computed from theprediction steps of the extended Kalman filter (EKF), unscented Kalmanfilter (UKF), Gauss von Mises (GVM) filter, and a particle filter andused in the accompanying EXAMPLE for demonstration according to oneembodiment of the present invention;

FIG. 10 is a graph showing plots of the normalized L₂ error usingdifferent methods for uncertainty propagation used in the accompanyingEXAMPLE for demonstration according to one embodiment of the presentinvention;

FIG. 11 is a block diagram illustrating components of an exemplaryoperating environment in which various embodiments of the presentinvention may be implemented; and

FIG. 12 is a block diagram illustrating an exemplary computer system inwhich embodiments of the present invention may be implemented.

DETAILED DESCRIPTION OF THE INVENTION

In the following description, for the purposes of explanation, numerousspecific details are set forth in order to provide a thoroughunderstanding of various embodiments of the present invention. It willbe apparent, however, to one skilled in the art that embodiments of thepresent invention may be practiced without some of these specificdetails. In other instances, well-known structures and devices are shownin block diagram form.

The ensuing description provides exemplary embodiments only, and is notintended to limit the scope, applicability, or configuration of thedisclosure. Rather, the ensuing description of the exemplary embodimentswill provide those skilled in the art with an enabling description forimplementing an exemplary embodiment. It should be understood thatvarious changes may be made in the function and arrangement of elementswithout departing from the spirit and scope of the invention as setforth in the appended claims.

Specific details are given in the following description to provide athorough understanding of the embodiments. However, it will beunderstood by one of ordinary skill in the art that the embodiments maybe practiced without these specific details. For example, circuits,systems, networks, processes, and other components may be shown ascomponents in block diagram form in order not to obscure the embodimentsin unnecessary detail. In other instances, well-known circuits,processes, algorithms, structures, and techniques may be shown withoutunnecessary detail in order to avoid obscuring the embodiments.

Also, it is noted that individual embodiments may be described as aprocess which is depicted as a flowchart, a flow diagram, a data flowdiagram, a structure diagram, or a block diagram. Although a flowchartmay describe the operations as a sequential process, many of theoperations can be performed in parallel or concurrently. In addition,the order of the operations may be re-arranged. A process is terminatedwhen its operations are completed, but could have additional steps notincluded in a figure. A process may correspond to a method, a function,a procedure, a subroutine, a subprogram, etc. When a process correspondsto a function, its termination can correspond to a return of thefunction to the calling function or the main function.

The terms “machine-readable medium” and/or “computer-readable medium”include, but are not limited to portable or fixed storage devices,optical storage devices, wireless channels and various other mediumscapable of storing, containing or carrying instruction(s) and/or data. Acode segment or machine-executable instructions may represent aprocedure, a function, a subprogram, a program, a routine, a subroutine,a module, a software package, a class, or any combination ofinstructions, data structures, or program statements. A code segment maybe coupled to another code segment or a hardware circuit by passingand/or receiving information, data, arguments, parameters, or memorycontents. Information, arguments, parameters, data, etc. may be passed,forwarded, or transmitted via any suitable means including memorysharing, message passing, token passing, network transmission, etc.

Furthermore, embodiments may be implemented by hardware, software,firmware, middleware, microcode, hardware description languages, or anycombination thereof. When implemented in software, firmware, middlewareor microcode, the program code or code segments to perform the necessarytasks may be stored in a machine readable medium. A processor(s) mayperform the necessary tasks.

Embodiments of the invention provide systems and methods for properlycharacterizing the uncertainty of a state defined on a n+1 dimensionalcylindrical manifold

^(n)×

, e.g., characterizing the uncertainty of the orbital state of anEarth-orbiting space object. More specifically, embodiments of thepresent invention provide a new class of multivariate probabilitydensity functions (PDFs), called Gauss von Mises (GVM) distributions,formulated to more accurately characterize the uncertainty of a spaceobject's state or orbit. Using the GVM distribution as input, extensionsand improvements can be made to tracking algorithms including theBayesian non-linear filter used for uncertainty propagation and datafusion, batch processing and orbit determination, multi-dimensionalquadrature, and likelihood ratios and other scoring metrics used in dataassociation. An exemplary application of embodiments described herein isthe use of the GVM distribution within a sequential non-linear filterresulting in improved algorithms for uncertainty propagation and datafusion to support higher level SSA functions. Embodiments of the presentinvention provide a tractable implementation of the Bayesian non-linearfilter using the GVM distribution to model state or orbital uncertaintywhile providing significant improvements, both in terms of accuracy andcomputational expense, over the aforementioned methods.

One feature of the GVM distribution is its definition on a cylindricalstate space

^(n)×

with the proper treatment of the angular coordinate within the generalframework of directional statistics. Hence, the GVM distribution isrobust for uncertainty quantification in orbital element space. TheGauss von Mises (GVM) distribution uses the von Mises distribution, theanalogy of a Gaussian distribution defined on a circle, to robustlydescribe uncertainty in the angular coordinate. The marginaldistribution in the Cartesian coordinates is Gaussian. Additionally, theGVM distribution can contain a parameter set controlling the correlationbetween the angular and Cartesian variables as well as the higher-ordercumulants which gives the level sets of the GVM PDF a distinctive“banana” or “boomerang” shape.

By providing a statistically robust treatment of the uncertainty in aspace object's orbital element state by rigorously defining theuncertainty on a cylindrical manifold, the GVM distribution can supporta suite of next-generation algorithms for uncertainty propagation, dataassociation, space catalog maintenance, and other SSA functions. Whenadapted to state PDFs modeled by GVM distributions, the general Bayesiannon-linear filter can be tractable. The prediction step of the resultingGVM filter can be made derivative-free, like the UKF, by new quadraturerules for integrating a function multiplied by a GVM weight function,thereby extending the unscented transform. Moreover, prediction usingthe GVM filter uses the propagation of the same number of sigma points(quadrature nodes) as the standard UKF. Thus, the GVM filter predictionstep (uncertainty propagation) has similar computational costs as theUKF and, as demonstrated in the EXAMPLE below, the former can maintain aproper characterization of the uncertainty much longer than the latter.In the most exceptional cases when the actual state uncertainty deviatesfrom a GVM distribution, a mixture version of the GVM filter can beformulated (using GVM distributions as the mixture components) toprovide proper uncertainty realism in analogy to the Gaussian sum(mixture) filter. Furthermore, the GVM filter can be applicable invarious regimes of space and to both conservative (e.g., gravity) andnon-conservative forces (e.g., atmospheric drag, solar radiationpressure). Stochastic process noise, uncertain model parameters, andresidual biases can also be treated within the GVM filter usingextensions of classical consider analysis and the Schmidt-Kalman filter.A maximum a posteriori batch processing capability for orbitdetermination (track initiation) can also be formulated which generatesa GVM PDF characterizing the initial orbital state and uncertainty froma sequence of input reports such as radar, electro-optical, or infraredsensor observation data or even full track states. To support the datafusion problem of tracking, the correction step of the Bayesiannon-linear filter can also be specialized to GVM distributions therebyenabling one to combine reports emanating from a common object toimprove the state or understanding of that object. The filter correctionstep can also furnish a statistically rigorous prediction error whichappears in the likelihood ratios for scoring the association of onereport to another. Thus, the new GVM filter can be used to supportmulti-target tracking within a general multiple hypothesis trackingframework. Additionally, the GVM distribution admits a distance metricwhich extends the classical Mahalanobis distance (χ² statistic). Thisnew “Mahalanobis von Mises” metric provides a test for statisticalsignificance and facilitates validation of the GVM filter. Anothernoteworthy feature of the GVM framework is its backwards compatibilitywith the space catalog and existing covariance-based algorithms. Thatis, the GVM distribution reduces to a Gaussian under suitable limits andthe GVM filter reduces to the UKF in the case of linear dynamical andmeasurement models.

To briefly demonstrate the implication of mis-characterizing orbitaluncertainty and the improvements obtained when using the GVMdistribution, FIGS. 9 a-9 f provides a simple example comparing theuncertainty propagation algorithms implicit in the EKF, UKF, and new GVMfilter. FIG. 9 a shows a particle representation of an orbital state PDFin equinoctial orbital element space plotted on the plane of thesemi-major axis and mean longitude coordinates. The particles 905 aredispersed according to the level curves 910 ranging from 0.5 to 3 sigmasin half sigma increments. FIG. 9 f shows the results of propagating thisinitial uncertainty (for eight orbital periods in this example) usingthe EKF, UKF, and GVM filter. The respective level curves of the PDFcomputed from the EKF 925, UKF 920, and GVM 915 filter are superimposedon the figure. The level curves computed from the GVM filter can capturethe actual uncertainty characterized by the particle ensemble. For theUKF, its covariance is indeed consistent (realistic) in the sense thatit agrees with that computed from the definition of the covariance usingthe true PDF. Thus, in this example, the UKF provides “covariancerealism” but does not support “uncertainty realism” since the covariancedoes not represent the actual banana-shaped uncertainty of the true PDF.Further, the state estimate produced from the UKF coincides with themean of the true PDF. However, the mean is displaced from the mode ofthe true PDF. Consequently, the probability that the object is within asmall neighborhood centered at the UKF state estimate (mean) isessentially zero. The EKF, on the other hand, provides a state estimatecoinciding closely with the mode of the true PDF, but the covariancetends to collapse making inflation necessary to begin to cover theuncertainty. In neither the EKF nor UKF case does the covarianceactually model the uncertainty. This example illustrates the problem ofusing a single Gaussian PDF (i.e., covariance) to represent uncertaintyand suggests that the GVM distribution provides a better representationof the actual orbital state PDF which is important if one is to achievestatistically robust characterization of uncertainty, which again isfundamental to achieving a robust capability across the SSN.

The Gauss von Mises (GVM) distribution and its resulting applicationscan provide a new suite of algorithms to support the future needs of theSSA mission. The distribution can provide a mechanism for moreaccurately characterizing the uncertainty in a space object's orbit atlittle or no additional cost to existing methods in operation whileproviding backwards compatibility with legacy systems. The framework cansupport (i) orbit determination to generate a realistic characterizationof an initial orbital state and uncertainty from a sequence of sensorobservations, (ii) non-linear uncertainty propagation to predict thefuture location of a space object and properly characterize its orbitaluncertainty at future times, (iii) conjunction analysis to estimate theprobability that two space objects will collide, (iv) anomaly andmaneuver detection to determine if an object deviates from its expectedkinematic trajectory, and (v) data fusion, tracking, and space catalogmaintenance to update a catalog of space objects, online, as newobservations are received, in order to improve the accuracy of theobjects in the catalog. These SSA tracking algorithms based on the GVMdistribution can provide next-generation upgrades to the AstrodynamicStandards software suite maintained by AFSpC and support the overallobjectives of the JSpOC Mission System. Various additional details ofembodiments of the present invention will be described below withreference to the figures.

The organization of the invention description is described below.

Section I entitled “Uncertainty on Manifolds” gives an overview ofuncertainty characterization in tracking and discusses coordinatesystems used to describe a space object's orbital state and the pitfallsof mistreating an angular coordinate by an unbounded Cartesiancoordinate. The von Mises distribution is then introduced as onepossible distribution to rigorously treat the uncertainty of an angularvariable.

Section II entitled “Construction of the Gauss von Mises Distribution”motivates and defines the Gauss von Mises (GVM) distribution.

Section III entitled “Elementary Properties” lists various mathematicaland statistical properties of the GVM distribution.

Section IV entitled “Gauss von Mises Quadrature” extends the methodologyof classical Gauss-Hermite quadrature to enable the computation of theexpected value of a non-linear transformation of a GVM random variable,thereby providing a general framework for GVM quadrature. Suchquadrature formulas play a key role in the derivation of the GVM filterprediction step.

Section V entitled “Uncertainty Propagation” develops the predictionstep of the Bayesian non-linear filter using the GVM distribution asinput thereby providing a method for approximating the non-lineartransformation of a GVM distribution as another GVM distribution.Performance metrics used for validation and extensions to thepropagation of mixtures of GVM distributions are also discussed.

Section VI entitled “Data Fusion” specializes the correction step of theBayesian non-linear filter to the case when the prior state is a randomvector represented by a GVM distribution and the update report is eitheranother GVM random vector or an observation related to the prior by astochastic measurement model. The corresponding prediction error termused to score the association of one report to another is also derived.

Section VII entitled “Batch Processing” treats the batch filteringproblem which processes a sequence of reports simultaneously to producea GVM distribution of the system state conditioned on the reports.

The reader of this document is assumed to be familiar with the standardnotation of linear algebra, in particular, the notions of matrices,vectors, matrix vector product, matrix inverse, and systems of linearequations. The reader of this document is also assumed to be familiarwith elementary probability theory, in particular, the notations ofrandom variables, random vectors, probability density functions, andtransformations of random variables and random vectors. The requiredbackground can be obtained by reading books associated with collegecourses in linear algebra, matrix analysis, and probability theory. Itmay also be useful to have familiarity with orbital mechanics.

I. Uncertainty on Manifolds

A permeating theme throughout space situational awareness is theachievement of the correct characterization and management ofuncertainty, which in turn is used to support conjunction analysis, dataassociation, anomaly detection, and sensor resource management. Thesuccess in achieving a proper characterization of the uncertainty in thestate of a space object can depend greatly on the choice of coordinatesystem. Under Gaussian assumptions, the coordinates used to representthe state space can impact how long one can propagate the uncertaintyunder a non-linear dynamical system. A Gaussian random vector need notget mapped to a Gaussian under a non-linear transformation. Therepresentation of a space object's kinematic state in orbital elementcoordinates, rather than Cartesian Earth-Centered-Inertial (ECI)position-velocity coordinates, is well-suited to the space surveillancetracking problem since such coordinates “absorb” the most dominant termin the non-linear gravitational force (i.e., the 1/r² term) leading to“more linear” propagations. Thus, these special coordinates can mitigatethe departure from “Gaussianity” under the non-linear propagation of aninitial Gaussian state probability density function (PDF) with respectto orbital elements. Additional discussions are provided in SubsectionI.A below. The benefits of the orbital element coordinates are alsoexploited in the propagation of a Gauss von Mises (GVM) PDF under thenon-linear two-body dynamics. This framework is developed later inSection V.

The application of orbital element coordinates within traditionalsequential filtering methods such as the extended Kalman filter (EKF),unscented Kalman filter (UKF), and even the highest fidelity Gaussiansum filters has one major disadvantage: the mean anomaly (or meanlongitude) angular coordinate describing the location along the orbit isincorrectly treated as an unbounded Cartesian coordinate. Some sideeffects and pitfalls of such mistreatments within the problems ofaveraging and fusing angular quantities are described in Subsection I.B.Ultimately, what is helpful to rectify these shortcomings is astatistically rigorous treatment of the uncertainty on the underlyingmanifold on which the system state is defined. The theory of directionalstatistics provides one possible development path. Though the theory cantreat uncertainty on very general manifolds (such as tori andhyper-spheres) possessing multiple directional quantities, this workfocuses on distributions defined on the circle

and the n+1-dimensional cylinder

^(n)×

. The latter is the manifold on which the orbital element coordinatesare more accurately defined. As a starting point, the von Mises PDF isintroduced in Subsection I.C as one possible distribution to rigorouslytreat the uncertainty of a single angular variable defined on thecircle. The von Mises distribution paves the way for the next sectionconcerning the development of the GVM distribution used to representuncertainty on a cylinder.

I.A. Orbital Element Coordinate Systems

With respect to Cartesian ECI position-velocity coordinates (r,{dot over(r)}), the acceleration {umlaut over (r)} of a space object (e.g.,satellite, debris) can be written in the form

$\begin{matrix}{\overset{¨}{r} = {{{- \frac{\mu_{\oplus}}{r^{3}}}r} + {{a_{pert}\left( {r,\overset{.}{r},t} \right)}.}}} & (1)\end{matrix}$

In this equation, r=|r|, μ_(⊕)=GM_(⊕) where G is the gravitationalconstant and M_(⊕) is the mass of the Earth, and a_(pert) encapsulatesall perturbing accelerations of the space object other than those due tothe two-body point mass gravitational acceleration.

The equinoctial orbital elements (a, h, k, p, q, l) define a system ofcurvilinear coordinates with respect to six-dimensionalposition-velocity space. Physical and geometric interpretations of thesecoordinates as well as the transformation from equinoctial elements toECI are known. Models for the perturbing acceleration a_(pert) are alsowell known.

The representation of the dynamical model (1) in coordinate systemsother than ECI can be obtained using the chain rule. Indeed, ifu=u(r,{dot over (r)}) denotes a coordinate transformation from ECIposition-velocity coordinates (r,{dot over (r)}) to a coordinate systemuε

⁶ (e.g., equinoctial orbital elements), then (1) is transformed to

$\begin{matrix}{{\overset{.}{u} = {{\overset{.}{u}}_{unpert} + {\frac{\partial u}{\partial\overset{.}{r}} \cdot {a_{pert}\left( {r,\overset{.}{r},t} \right)}}}},{where}} & (2) \\{{\overset{.}{u}}_{unpert} = {{\frac{\partial u}{\partial r} \cdot \overset{.}{r}} - {\frac{\mu_{\oplus}}{r^{3}}{\frac{\partial u}{\partial\overset{.}{r}} \cdot {r.}}}}} & (3)\end{matrix}$

If u is the vector of equinoctial orbital elements, then (3) simplifiesto

{dot over (u)} _(unpert)=(0,0,0,0,0,√{square root over(μ_(⊕)/α³)})^(T).  (4)

Therefore, the time evolution of a space object's equinoctial orbitalelement state (a, h, k, p, q, l) under the assumption of unperturbedtwo-body dynamics is

a(t)=a ₀ ,h(t)=h ₀ ,k(t)=k ₀ ,p(t)=p ₀ ,q(t)=q ₀ ,l(t)=l ₀ +n ₀(t−t₀),  (5)

where n₀=√{square root over (μ_(⊕)/a₀ ³)} is the mean motion at theinitial epoch.

It is now argued that an equinoctial orbital element state can beregarded as a state defined on the six-dimensional cylinder

⁵×

. The definition of the equinoctial orbital elements (a, h, k, p, q, l)implicitly assumes closed (e.g., elliptical) orbits which imposes thefollowing constraints:

a>0,h ² +k ²<1,p,qε

,lε

.

Though any real value can be assigned to the mean longitude coordinatel, the angles l and l+2πk, for any integer kε

, define the same location along the orbit. Hence, l can be regarded asan angular coordinate defined on the circle

. Mathematically, the first five elements (a, h, k, p, q) are notdefined everywhere on

⁵. Physically however, it can be argued that the sample space in thesefive coordinates can be extended to all of

⁵. Indeed, the elements (a, h, k, p, q) describe the geometry of the(elliptical) orbit and the orientation of the orbit relative to theequatorial plane. In practice, the uncertainties in the orbital geometryand orientation are sufficiently small so that the probability that anelement is close to the constraint boundary is negligible. Moreover,under unperturbed dynamics, these five elements are conserved byKepler's laws (see Equations (5)) and, consequently, their respectiveuncertainties do not grow. That said, under the above assumptions, themanifold on which the elements (a, h, k, p, q) are defined can beapproximated as all of

⁵ and the full six-dimensional equinoctial orbital element state spacecan be defined on the cylinder

⁵×

. With the inclusion of perturbations in the dynamics (e.g.,higher-order gravity terms, drag, solar radiation, etc.), the first fiveelements evolve with small periodic variations in time and theiruncertainties exhibit no long-term secular growth. Thus, the samecylindrical assumptions on the state space apply. These discussions alsoequally apply to other systems of orbital elements such as Poincaréorbital elements, modified equinoctial orbital elements, and alternateequinoctial orbital elements.

Though the uncertainties in the first five equinoctial elements exhibitonly small periodic changes under two-body dynamics, it is theuncertainty along the semi-major axis coordinate a which causes theuncertainty along the mean longitude coordinate l to grow without bound.In other words, as time progresses, one can generally maintain a goodunderstanding of the geometry and orientation of the orbit, butconfidence is gradually lost in the exact location of the object alongits orbit. The growth in the uncertainty in l can cause undesirableconsequences if l is incorrectly treated as an unbounded Cartesianvariable or the uncertainty in l is modeled as a Gaussian. With aGaussian assumption imposed on the PDF in l, one would assign differentlikelihoods to l and l+2πk for different values of kε

, even though the mean longitudes l and l+2πk define the same locationon the circle. Additional insights are provided in the next subsection.

I.B. Examples of Improper Treatments of Angular Quantities

As described in the previous subsection, the mean longitude coordinate(i.e., the sixth equinoctial orbital element) is an angular coordinatein which the angles l and l+2πk are identified (equivalent) for anyinteger kε

. In other words, the mean longitude is a circular variable and, assuch, a rigorous treatment should not treat it as an unboundedreal-valued Cartesian variable. In many cases (and with proper branchcuts defined), it is practical to treat the mean anomaly as real-valued,allowing for distributions to be represented as Gaussians, for example.The drawback of this approach is that the resulting statistics candepend on choice of the integer ‘k’.

For example, suppose one wishes to compute the conventional average θ oftwo angles θ₁=π/6 and θ₂=−π/6. On one hand,

$\overset{\_}{\theta} = {{\frac{1}{2}\left( {\theta_{1} + \theta_{2}} \right)} = 0.}$

On the other hand, since θ₂=−π/6 and θ₂=11π/6 define the same locationon the circle, then using this equivalent value of θ₂ in the definitionof the average would yield θ=π. Thus, in order for the conventionalaverage to be well-defined, the statistic

$\overset{\_}{\theta} = {{\frac{1}{2}\left\lbrack {\left( {\theta_{1} + {2\pi \; k_{1}}} \right) + \left( {\theta_{2} + {2\pi \; k_{2}}} \right)} \right\rbrack} = {{\frac{1}{2}\left( {\theta_{1} + \theta_{2}} \right)} + {\pi \left( {k_{1} + k_{2}} \right)}}}$

must be independent of the choice of k₁, k₂ε

. Clearly this is not the case for the simple example studied above.Different values of k₁ and k₂ can define averages at opposite sides ofthe circle. The unscented transform used in the UKF can also exhibitsimilar side effects when used in equinoctial orbital element space(especially if the mean longitude components of the sigma points aresufficiently dispersed) because it requires one to compute a weightedaverage of angular (mean longitude) components.

A more striking example of the mis-representation of a circular variableas an unbounded Cartesian variable is demonstrated in FIG. 6 a. In thissimple example, consider a one-dimensional angular state space with twoindependent states: a “prior” θ₁ and an “update” θ₂ with respective PDFsp₁(θ₁) and p₂(θ₂). The PDF in θ₁ is diffuse (i.e., the uncertainty islarge) and, further, θ₁ is incorrectly modeled as an unbounded Cartesianvariable. The update state θ₂ has a mean of zero and a very smallvariance (uncertainty). Any one of these curves can represent the PDF ofthe update since they are all equivalent up to an integer 2π shift. The“correct” update is ambiguous. Within a non-linear filteringapplication, one might want to fuse the prior θ₁ with the informationfrom θ₂. In such a case, the Bayesian filter correction step yields a“fused PDF” p_(ƒ)(θ) in the fused state θ given by

$\begin{matrix}{{{p_{f}(\theta)} = {\frac{1}{c}{p_{1}(\theta)}{p_{2}(\theta)}}},{c = {\int_{I}^{\;}{{p_{1}(\theta)}{p_{2}(\theta)}{\theta}}}},} & (6)\end{matrix}$

where I=(−∞, ∞). Thus, the expression obtained for the fused PDF as wellas the normalization constant c depend on the choice (i.e., 2π shift) ofthe update PDF. In particular, a mis-computed normalization constant ccan have severe consequences within a tracking system, as analogousexpressions for c appear in the likelihood ratios for scoring theassociation of one report to another. FIG. 6 b depicts how theseambiguity issues can be resolved by turning to the theory of directionalstatistics. Further analysis is provided in the next subsection.

I.C. von Mises Probability Density Function

The source of the ambiguity in computing the fused PDF (6) in theexample of the previous subsection lies in the incorrect treatment of anangular variable as an unbounded Cartesian real-valued variable. Theseproblems can be rectified by representing the state PDF not necessarilyon a Cartesian space

^(n) but instead on the manifold that more accurately describes theglobal topology of the underlying state space. One change made for(equinoctial) orbital element states can be the representation of theuncertainty in the angular coordinate (mean longitude) as a PDF on thecircle

so that the joint PDF in the six orbital elements defines a distributionon the cylinder

⁵×

. In order to do so, PDFs defined on the circle are used.

A function p:

→

is a probability density function (PDF) on the circle

if

1. p(θ)≧0 almost everywhere on (−∞, ∞),

2. p(θ+2π)=p(θ) almost everywhere on (−∞, ∞),

3. ∫_(−π) ^(π)p(θ)dθ=1.

One example satisfying the above conditions is the von Misesdistribution which provides an analogy of a Gaussian distributiondefined on a circle. The von Mises PDF in the angular variable θ isdefined by

$\begin{matrix}{{{\left( {{\theta;\alpha},\kappa} \right)} = \frac{^{\kappa \; {\cos {({\theta - \alpha})}}}}{2\pi \; {I_{0}(\kappa)}}},} & (7)\end{matrix}$

where I₀ is the modified Bessel function of the first kind of order 0.The parameters μ and κ are measures of location and concentration. Asκ→0⁺, the von Mises distribution tends to a uniform distribution. Forlarge κ, the von Mises distribution becomes concentrated about the angleα and approaches a Gaussian distribution in θ with mean μ and variance1/κ. Some example plots are shown in FIG. 7. An algebraically equivalentexpression of (7) which is numerically stable for large values of κ is

$\begin{matrix}{{\left( {{\theta;\alpha},\kappa} \right)} = {\frac{^{{- 2}\kappa \; \sin^{2}\frac{1}{2}{({\theta - \alpha})}}}{2\pi \; ^{- \kappa}{I_{0}(\kappa)}}.}} & (8)\end{matrix}$

The von Mises distribution (8) can be used in the next section toconstruct the Gauss von Mises distribution defined on the cylinder

^(n)×

.

This section concludes with a discussion on how the von Misesdistribution can resolve the problems in the previous subsectionconcerning the averaging and fusion of angular quantities. Suppose θ₁, .. . , θ_(N) is a sample of independent observations coming from a vonMises distribution. Then, the maximum likelihood estimate of thelocation parameter α is

$\begin{matrix}{\hat{\alpha} = {{\arg \left( {\frac{1}{N}{\sum\limits_{j = 1}^{N}^{i\; \theta_{j}}}} \right)}.}} & (9)\end{matrix}$

Unlike the conventional average (θ₁+ . . . +θ_(N))/N, the “average”{circumflex over (α)} is independent of any 2π shift in the anglesθ_(j). In the fusion example, if the prior and update PDFs are von Misesdistributions, i.e.,

p ₁(θ₁)=

(θ₁;α₁,κ₁),p ₂(θ₂)=

(θ₂;α₂,κ₂),

then the fused PDF (6) can be formed unambiguously andp_(ƒ)(θ)=p_(ƒ)(θ+2πk) for any kε

. As seen in FIG. 6 b, there is no longer any ambiguity in how to choosethe “correct” update PDF p₂(θ), since both p₁(θ) and p₂(θ) are properlydefined circular PDFs. The integration interval I in the computation ofthe normalization constant c in (6) is any interval of length 2π, i.e.,I=(a,a+2π) for some aε

. In other words, the value of c is independent of the choice of a.

II. Construction of the Gauss von Mises Distribution

This section defines the Gauss von Mises (GVM) distribution used tocharacterize the uncertainty in a space object's orbital state. Later,in Subsection II.B, the connection between the GVM distribution and arelated family of distributions is discussed.

The proposed GVM distribution is not defined as a function of otherrandom variables with specified probability density functions (PDFs).Instead, its construction is based on satisfying the followingconditions listed below.

-   -   1. The GVM family of multivariate PDFs can be defined on the        n+1-dimensional cylindrical manifold        ^(n)×        .    -   2. The GVM distribution can have a sufficiently general        parameter set which can model non-zero higher-order cumulants        beyond a usual “mean” and “covariance.”    -   3. The GVM distribution can account for correlation between the        Cartesian random vector xε        ^(n) and the circular random variable θε        .    -   4. The level sets of the GVM distribution can be generally        “banana” or “boomerang” shaped so that a more accurate        characterization of the uncertainty in a space object's state in        equinoctial orbital elements can be achieved.    -   5. The GVM distribution can reduce to a multivariate Gaussian        PDF in a suitable limit.    -   6. The GVM distribution can permit a tractable implementation of        the general Bayesian nonlinear filter and other applications to        support advanced SSA.

In clarification of the first condition, a function p:

^(n)×

→

is a probability density function (PDF) on the cylinder

^(n)×

if

1. p(x,θ)>0 almost everywhere on

_(n)×

,

2. p(x,θ+2π)=p(x,θ) almost everywhere on

^(n)×

,

3.

_(n)∫_(−π) ^(π)p(x,θ)dθdx=1.

This definition extends the definition of a PDF defined on a circle. Oneexample satisfying the above conditions is

p(x,θ)=

(x;μ,P)

(θ;α,κ),

where

(θ;α,κ) is the von Mises PDF defined by (8) and

(x;μ,P) is the multivariate Gaussian PDF given by

$\begin{matrix}{{{\left( {{x;\mu},P} \right)} = {\frac{1}{\sqrt{\det \left( {2\pi \; P} \right)}}{\exp \left\lbrack {{- \frac{1}{2}}\left( {x - \mu} \right)^{T}{P^{- 1}\left( {x - \mu} \right)}} \right\rbrack}}},} & (10)\end{matrix}$

where με

^(n) and P is an n×n symmetric positive-definite (covariance) matrix.This simple example, though used as a starting point in constructing theGVM distribution, does not satisfy all of the conditions listed earlier.In particular, the desired family of multivariate PDFs should modelcorrelation between x and θ and have level sets possessing a distinctivebanana or boomerang shape.

To motivate the construction of the GVM distribution, consider tworandom vectors xε

^(n) and yε

^(m) whose joint PDF is Gaussian:

$\begin{matrix}{{p\left( {x,y} \right)} = {{\left( {{{\begin{bmatrix}x \\y\end{bmatrix};}\begin{bmatrix}\mu_{x} \\\mu_{y}\end{bmatrix}},\begin{bmatrix}P_{xx} & P_{yx}^{T} \\P_{yx} & P_{yy}\end{bmatrix}} \right).}}} & (11)\end{matrix}$

Using the definition of conditional probability and the Schur complementdecomposition, the joint PDF (11) can be expressed as

$\begin{matrix}{\begin{matrix}{{p\left( {x,y} \right)} = {{p(x)}{p\left( y \middle| x \right)}}} \\{= {\left( {{x;\mu_{x}},P_{xx}} \right)\begin{pmatrix}{{y;{\mu_{y} + {P_{yx}{P_{xx}^{- 1}\left( {x - \mu_{x}} \right)}}}},} \\{P_{yy} - {P_{yx}P_{xx}^{- 1}P_{yx}^{T}}}\end{pmatrix}}} \\{{\equiv {\left( {{x;\mu_{x}},P_{xx}} \right)\left( {{y;{g(x)}},Q} \right)}},}\end{matrix}{where}{{{g(x)} = {\mu_{y} + {P_{yx}{P_{xx}^{- 1}\left( {x - \mu_{x}} \right)}}}},{Q = {P_{yy} - {P_{yx}P_{xx}^{- 1}{P_{yx}^{- 1}.}}}}}} & (12)\end{matrix}$

One very powerful observation resulting from this Schur complementdecomposition is that (12) defines a PDF in (x,y) for an (analytic)function g(x) and a symmetric positive-definite matrix Q. In particular,if y is univariate and relabelled as θ, then

p(x,θ)=

(x;μ _(x) ,P _(xx))

(θ;Θ(x),κ),

defines a PDF on

^(n)×

for an analytic function Θ:

^(n)→

and a positive scalar κ. To make it robust for a circular variable θε

and hence define a PDF on the cylinder

^(n)×

, the Gaussian PDF in θ can be replaced by the von Mises PDF in θ:

p(x,θ)=

(x;μ _(x) ,P _(xx))

(θ;Θ(x),κ).  (13)

The definition of the Gauss von Mises distribution fixes the specificform of the function Θ(x) in (13) so that the PDF can model non-zerohigher-order cumulants (i.e., the banana shape of the level sets) but isnot overly complicated so as the make the resulting Bayesian filterprediction and correction steps intractable.

Definition [Gauss von Mises (GVM) Distribution] The random variables(x,θ)ε

^(n)×

are said to be jointly distributed as a Gauss von Mises (GVM)distribution if their joint probability density function has the form

$\begin{matrix}{{p\left( {x,\theta} \right)} = {\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right)}} \\{{\equiv {\left( {{x;\mu},P} \right)\left( {{\theta;{\Theta (x)}},\kappa} \right)}},}\end{matrix}$ where${{\left( {x,{;\mu},P} \right)} = {\frac{1}{\sqrt{\det \left( {2\pi \; P} \right)}}{\exp \left\lbrack {{- \frac{1}{2}}\left( {x - \mu} \right)^{T}{P^{- 1}\left( {x - \mu} \right)}} \right\rbrack}}},{{\left( {{\theta;{\Theta (x)}},\kappa} \right)} = {\frac{1}{2\pi \; e^{- \kappa}{I_{0}(\kappa)}}{\exp \left\lbrack {{- 2}\kappa \; \sin^{2}\frac{1}{2}\left( {\theta - {\Theta (x)}} \right)} \right\rbrack}}},{{{and}{\Theta (x)}} = {\alpha + {\beta^{T}z} + {\frac{1}{2}z^{T}\Gamma \; z}}},{z = {A^{- 1}\left( {x - \mu} \right)}},{P = {{AA}^{T}.}}$

The parameter set (μ, P, α, β, Γ, κ) can be subject to the followingconstraints: με

^(n), P can be an n×n symmetric positive-definite matrix, αε

, βε

^(n), Γ can be an n×n symmetric matrix, and κ≧0. The matrix A in thedefinition of the normalized variable z can be the lower-triangularCholesky factor of the parameter matrix P. The notation

(x,θ)˜GVM(μ,P,α,β,Γ,κ)

denotes that (x,θ) can be jointly distributed as a GVM distribution withthe specified parameter set.

It is noted that the function Θ(x) appearing in the GVM distribution canbe an inhomogeneous quadratic in x or, equivalently, the normalizedvariable z. The use of the normalized variable z in the definition canbe made to simplify the ensuing mathematics and so that the parameters βand Γ are dimensionless. Additionally, the parameters β and Γ modelcorrelation between x and θ. The parameter matrix Γ can be tuned to givethe level sets of the GVM distribution their distinctive banana orboomerang shape. In subsequent sections, it is demonstrated that thedefinition of the GVM distribution satisfies the remaining conditionslisted at the beginning of the section.

II.A. Relationship to Mardia and Sutton

A related distribution defined on a cylindrical manifold proposed byMardia and Sutton can be recovered starting from Equation (12). If x isunivariate and relabelled as θ, then

p(y,θ)=

(y;g(θ),Q)

(θ;μ_(θ) ,P _(θθ))

defines a PDF on

^(m)×

for any analytic function g:

→

^(m). A PDF defined on the cylinder

^(m)×

can be obtained by replacing the Gaussian PDF in θ by the von Mises PDFin θ (and relabelling the other parameters) to yield

p(y,θ)=

(y;g(θ),P)

(θ;α,κ).  (14)

One specific choice for the function g(θ) is

g(θ)=a ₀ +a ₁ cos θ+b ₁ sin θ,  (15)

for some parameter vectors a₀, a₁, b₁ε

^(m).

The PDF (14) with g(θ) specified by (15) can be a generalization of thebivariate distribution defined on the two-dimensional cylinder

×

proposed by Mardia and Sutton. Their derivation also differs from thatpresented here in the sense that the authors start with a trivariateGaussian distribution and then “project it” onto the cylinder.

Though the Mardia-Sutton distribution can be used in place of the GVMdistribution in estimation and other space situational awarenessalgorithms, it is noted that the parameter set in the former providesless control over the magnitude of the higher-order cumulants and thebending of the banana or boomerang shaped level sets. For these reasons,the GVM distribution is preferred since it overcomes these limitations.

III. Elementary Properties

This section lists mathematical and statistical properties of the Gaussvon Mises (GVM) distribution based on the definition stated in SectionII and provides an outline of the derivation of non-trivial properties.These elementary properties are applied in subsequent sections of thedescription.

III.A. Mode

The GVM distribution is unimodal on

^(n)×

with

$\left( {\mu,\alpha} \right) = {\underset{x,\theta}{argmax}{\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right).}}$

III.B. Characteristic Function

The characteristic function of the GVM distribution, i.e.,

ϕ  ( ξ , m ; μ , P , α , β , Γ , κ ) ≡  E  [  i  ( ξ T  x + m  θ ) ] =  ∫ ℝ n  ∫ - π π   i  ( ξ T  x + m   θ )   ( x , θ ; μ, P , α , β , Γ , κ )   θ   x ,   for   ξ ∈ ℝ n  m ∈ ℤ   is  ϕ  ( ξ , m ; μ , P , α , β , Γ , κ ) = 1 det  ( I - im   Γ )  I m   ( κ ) I 0  ( κ ) × exp  [ i  ( μ T  ξ + m   α ) - 1 2  ( AT  ξ + m   β ) T  ( I - im   Γ ) - 1  ( A T  ξ + m   β ) ] , (16 )

where I denotes the n×n identity matrix and I_(p) is the modified Besselfunction of the first kind of order p.

To derive (16), the definition of the GVM distribution in conjunctionwith the definition of the characteristic function is applied to obtain

ϕ  ( ξ , m ; μ , P , α , β , Γ , κ ) =  ∫ ℝ n   i   ξ T  x   (x ; μ , P )  [ ∫ - π π   im   θ   ( θ ; Θ  ( x ) , κ )   θ ]  x =  ∫ ℝ n   i   ξ T  x   ( x ; μ , P )  I  m   ( κ ) I 0 ( κ )   im   Θ  ( x )   x ,

This integral can be expressed in the form

ϕ  ( ξ , m ; μ , P , α , β , Γ , κ ) = 1 det  ( 2  π   P )  I  m  ( κ ) I 0  ( κ )  ∫ ℝ n   g  ( x )   x ,  where   g  ( x) =  - 1 2  ( x - μ ) T  P - 1  ( x - μ ) + im   Θ  ( x ) + i  ξ T  x ≡  g  ( x * ) - 1 2  ( x - x * ) T  A - T  ( I - im   Γ ) A - 1  ( x - x * ) ,   and   x * =  μ + iA  ( I - im   Γ ) -1  ( A T  ξ + m   β ) , g  ( x * ) =  i  ( μ T  ξ + m   α ) -1 2  ( A T  ξ + m   β ) T  ( I - im   Γ ) - 1  ( A T  ξ + m  β ) .

This decomposition of g(x) can follow by “completing the square.” Theresult (16) follows from known integration formulas for integrandscontaining an exponential of a quadratic form and elementary propertiesof the determinant.

III.C. Low-Order Moments

If (x,θ)˜GVM(μ, P, α, β, Γ, κ), then

${{E\lbrack x\rbrack} = \mu},{{E\left\lbrack {\left( {x - \mu} \right)\left( {x - \mu} \right)^{T}} \right\rbrack} = P},{{E\left\lbrack ^{\theta} \right\rbrack} = {\frac{1}{\sqrt{\det \left( {I - {\Gamma}} \right)}}\frac{I_{1}(\kappa)}{I_{0}(\kappa)}{\exp \left\lbrack {{\alpha} - {\frac{1}{2}{\beta^{T}\left( {I - {\Gamma}} \right)}^{- 1}\beta}} \right\rbrack}}},{{E\begin{bmatrix}^{\theta} & z\end{bmatrix}} = {{\left( {I - {\Gamma}} \right)}^{- 1}\beta \; {E\left\lbrack ^{\theta} \right\rbrack}}},{{E\begin{bmatrix}^{\theta} & {zz}^{T}\end{bmatrix}} = {\left\lbrack {\left( {I - {\Gamma}} \right)^{- 1} - {\left( {I - {\Gamma}} \right)^{- 1}{{\beta\beta}^{T}\left( {I - {\Gamma}} \right)}^{- 1}}} \right\rbrack {E\left\lbrack ^{\theta} \right\rbrack}}},$

where z=A⁻¹(x−μ) and P=AA^(T). These results follow from thecharacteristic function (16).

III.D. Differential Entropy

If (x,θ)˜GVM(μ, P, α, β, Γ, κ), then its differential entropy can be

$\begin{matrix}\begin{matrix}{{H\left( {x,\theta} \right)} = {E\left\lbrack {{- \ln}\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right)} \right\rbrack}} \\{= {{\frac{1}{2}\ln \mspace{14mu} {\det \left( {2\pi \; {eP}} \right)}} + {\ln \left( {2\pi \; {I_{0}(\kappa)}} \right)} - {\kappa \; {\frac{I_{1}(\kappa)}{I_{0}(\kappa)}.}}}}\end{matrix} & (17)\end{matrix}$

III.E. Invariance Property

If (x,θ)˜GVM(μ, P, α, β, Γ, κ) defined on

^(n)×

and

$\begin{matrix}{{\overset{\sim}{x} = {{Lx} + d}},{\overset{\sim}{\theta} = {\theta + \alpha + {b^{T}x} + {\frac{1}{2}x^{T}{Cx}}}},} & (18)\end{matrix}$

where L is an m×n matrix (with m≦n), dε

^(m), aε

, bε

^(n), and C is a symmetric n×n matrix, then

({tilde over (x)},{tilde over (θ)})˜GVM({tilde over (μ)},{tilde over(P)},{tilde over (α)},{tilde over (β)},{tilde over (Γ)},{tilde over(κ)})

defined on

^(m)×

, where

$\begin{matrix}{{\overset{\sim}{\mu} = {{L\; \mu} + d}},{\overset{\sim}{P} = {L\; P\; L^{T}}},{\overset{\sim}{\alpha} = {\alpha + a + {b^{T}\mu} + {\frac{1}{2}\mu^{T}C\; \mu}}},{\overset{\sim}{\beta} = {Q^{T}\left\lbrack {\beta + {A^{T}\left( {{C\; \mu} + b} \right)}} \right\rbrack}},{\overset{\sim}{\Gamma} = {{Q^{T}\left( {\Gamma + {A^{T}{CA}}} \right)}Q}},{\overset{\sim}{\kappa} = {\kappa.}}} & (19)\end{matrix}$

In these equations, Q can be an n×m matrix with mutually orthonormalcolumns and R can be an m×m positive-definite upper-triangular matrixsuch that QR=(LA)^(T). (The Q and R matrices can be computed byperforming a “thin” QR factorization of (LA)^(T).) Additionally, thelower-triangular Cholesky factor Ã such that {tilde over (P)}=ÃÃ^(T) isÃ=R^(T).

This result, which states that a GVM distribution remains a GVMdistribution under transformations of the form (18), follows uponcomputing the characteristic function of the transformed randomvariables ({tilde over (x)},{tilde over (θ)}) (in analogy to thederivation in Subsection III.B) and then identifying it with thecharacteristic function (16) possessing the parameters defined inEquations (19).

III.F. Transformation to Canonical Form

The transformation

$\begin{matrix}{{z = {A^{- 1}\left( {x - \mu} \right)}},{\varphi = {\theta - \alpha - {\beta^{T}z} - {\frac{1}{2}z^{T}\Gamma \; z}}},} & (20)\end{matrix}$

reduces (x,θ)˜GVM(μ, P, α, β, Γ, κ) to the canonical or standardized GVMdistribution with PDF

p(z,φ)=

(z,φ;0,I,0,0,0,κ)=

(z;0,I)

(φ;0,κ).

This result follows by noting that the transformation (20) is a specialcase of (18).

III.G. Marginalization

If (x,θ)˜GVM(μ, P, α, β, Γ, κ), then the marginal distribution of x isthe Gaussian PDF

(x; μ, P). This result follows immediately from the definition of theGVM distribution:

$\begin{matrix}{{p(x)} = {\int_{- \pi}^{\pi}{\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right){\theta}}}} \\{= {{\left( {{x;\mu},P} \right)}{\int_{\pi}^{\pi}{{{\mathcal{M}}\left( {{\theta;{\Theta (x)}},\kappa} \right)}{\theta}}}}} \\{= {{\left( {{x;\mu},P} \right).}}}\end{matrix}$

The marginal distribution of ({tilde over (x)},θ), where {tilde over(x)}=(x_(i) ₁ , . . . , x_(i) _(m) )^(T) and 1≦m<n=dim(x), is

p({tilde over (x)},θ)=

({tilde over (x)},θ;{tilde over (μ)},{tilde over (P)},{tilde over(α)},{tilde over (β)},{tilde over (Γ)},{tilde over (κ)}),

where the transformed parameter set is specified by Equations (19) witha=0, b=0, C=0, d=0, and

$L = {\begin{bmatrix}e_{i_{1}}^{T} \\\vdots \\e_{i_{m}}^{T}\end{bmatrix}.}$

In this construction of the matrix L, e_(j) is the n×1 vector where thej-th component is unity and all other components are zero.

III.H. Conversion between Gauss von Mises and Gaussian Distributions

As κ→∞ and ∥Γ∥→0, a GVM distribution in (x,θ) becomes a Gaussiandistribution in x=(x,θ). Under these conditions,

$\begin{matrix}{{{\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right)} \approx {\left( {{x;},} \right)}},{where}} & (21) \\{{ = \begin{bmatrix}\mu \\\alpha\end{bmatrix}},{= {}^{T}},{= {\begin{bmatrix}A & 0 \\\beta^{T} & {1/\sqrt{\kappa}}\end{bmatrix}.}}} & (22)\end{matrix}$

The approximation (21) is justified by first noting that, for fixed x,the von Mises distribution

(θ; Θ(x), κ) becomes concentrated at the point θ=Θ(x) as κ→∞. Thus, forlarge κ,

$\begin{matrix}{{\left( {{\theta;{\Theta (x)}},\kappa} \right)} = {\frac{1}{2{\pi }^{- \kappa}{I_{0}(\kappa)}}{\exp \left\lbrack {{- 2}\kappa \; \sin^{2}\frac{1}{2}\left( {\theta - {\Theta (x)}} \right)} \right\rbrack}}} \\{\approx {\sqrt{\frac{\kappa}{2\pi}}{\exp \left\lbrack {{- \frac{1}{2}}{\kappa \left( {\theta - {\Theta (x)}} \right)}^{2}} \right\rbrack}}} \\{{= {\left( {{\theta;{\Theta (x)}},{1/\kappa}} \right)}},}\end{matrix}$

noting that

${{{{\left. {I_{0}(\kappa)} \right.\sim ^{\kappa}}/\sqrt{2{\pi\kappa}}}\mspace{14mu} {as}\mspace{14mu} \kappa}->{{\infty \mspace{14mu} {and}\mspace{14mu} \sin^{2}\frac{1}{2}{\left. \varphi \right.\sim\frac{1}{4}}\varphi^{2}\mspace{14mu} {as}\mspace{14mu} \varphi}->0}};$

hence

(x,θ;μ,P,α,β,Γ,κ)≈

(x;μ,P)

(θ;Θ(x),1/κ).  (23)

As ∥Γ∥→0, Θ(x) tends to a linear function in x thereby reducing theright-hand side of (23) to a Gaussian in x and θ.

The expressions for the mean

and covariance

in (22) of the approximating Gaussian in

=(x,θ) are derived as follows. First, define

q

(x,θ)=−ln

(x,θ;μ,P,α,β,Γ,κ),

(

)=−ln

(

;

,

).

The approximating Gaussian in (21) is the “osculating Gaussian” definedas the Gaussian which is tangent to the GVM distribution at the mode.Specifically, the mean

is the mode of the GVM distribution, as specified in (22), and thecovariance

is obtained by demanding equality between the second-order partialderivatives of

and

evaluated at the mode. Indeed,

∂ 2  q  ∂ x   ∂ x T   x = m =  - 1 ,  ∂ 2  q ℳ ∂ x   ∂ x T  ( x , θ ) = ( μ , α ) = A - T  ( I + κββ T )  A - 1 ,  ∂ 2  q ∂θ   ∂ x   ( x , θ ) = ( μ , α ) = - κ   A - T  β ,  ∂ 2  q ℳ∂ θ 2   ( x , θ ) = ( μ , α ) = κ ;   hence   - 1 = [ A - T  (I + κββ T )  A - 1 - κ   A - T  β - κβ T  A - 1 κ ] . ( 24 )

Inverting (24) using the block matrix inversion lemma yields

$\begin{matrix}{= {\begin{bmatrix}{AA}^{T} & {A\; \beta} \\{\beta^{T}A^{T}} & {{\beta^{T}\beta} + {1/\kappa}}\end{bmatrix}.}} & (25)\end{matrix}$

Finally, computing the lower-triangular Cholesky factor

of the covariance (25) results in the expression in (22). It is notedthat the expressions in (22) can also be recovered from (23) by settingΓ=0 and combining the two Gaussian PDFs into a single Gaussian PDF in

by way of a “completing the square” operation.

In summary, Equations (21) and (22) indicate how one can convert to andfrom a GVM distribution and a Gaussian. The approximation of a Gaussiandistribution in x=(x,θ) by a GVM distribution with Γ=0 becomes exact inthe limit as the standard deviation in θ tends to zero. Theapproximation of a GVM distribution by a Gaussian becomes exact as κ→∞and ∥Γ∥→0. In any case, the equations provide an approximation of a GVMdistribution by the osculating or tangent Gaussian distribution at themodal point.

III.I Mahalanobis von Mises Statistic

Given (x,θ)˜GVM(μ, P, α, β, Γ, κ), the Mahalanobis von Mises statisticis defined by

$\begin{matrix}{{{M\left( {x,{\theta;\mu},P,\alpha,Β,\Gamma,\kappa} \right)} \equiv {{\left( {x - \mu} \right)^{T}{P^{- 1}\left( {x - \mu} \right)}} + {4\kappa \; \sin^{2}\frac{1}{2}\left( {\theta - {\Theta (x)}} \right)}}},} & (26)\end{matrix}$

where

${{\Theta (x)} = {\alpha + {\beta^{T}z} + {\frac{1}{2}z^{T}\Gamma \; z}}},{z = {A^{- 1}\left( {x - \mu} \right)}},{and}$P = AA^(T).

For large κ, M{dot over (˜)}χ²(n+1), where χ²(ν) is the chi-squaredistribution with v degrees of freedom and ‘{dot over (˜)}’ means“approximately distributed.” The statistic (26) is analogous to theMahalanobis distance for a multivariate Gaussian random vector x˜N(μ,P). Indeed, the first term in (26) is precisely this Mahalanobisdistance. Like the classical Mahalanobis distance, the Mahalanobis vonMises statistic provides a mechanism for testing if some realization(x_(*),θ_(*)) from a GVM distribution is statistically significant.Additional details on the interpretation of the Mahalanobis von Misesstatistic are provided following the derivation below.

The derivation of the distribution of M given that (x,θ) are jointlydistributed as a GVM distribution is as follows. Applying thetransformation (20) to canonical form yields

${M = {{z^{T}z} + {4\kappa \; \sin^{2}\frac{1}{2}\varphi}}},$

where now z˜N(0,I) and φ˜VM(0,κ) with z and φ independent. Therefore,

M=χ ²(n)+Y _(κ),

where

$Y_{\kappa} \equiv {4\; \kappa \; \sin^{2}\frac{1}{2}\varphi}$

and is independent of the χ²(n) random variable. It follows fromelementary probability theory that the PDF of Y_(κ) is

$\begin{matrix}{{p_{Y_{\kappa}}(y)} = \left\{ \begin{matrix}{\frac{^{{- y}/2}}{\pi \; ^{- \kappa}{I_{0}(\kappa)}\sqrt{y\left( {{4\; \kappa} - y} \right)}},} & {{0 < y < {4\; \kappa}},} \\{0,} & {{otherwise}.}\end{matrix} \right.} & (27)\end{matrix}$

Noting that I₀(κ)˜r^(κ)/√{square root over (2πκ)} and √{square root over(y(4κ−y))}˜2√{square root over (κy)} as κ→∞, it follows that

${\left. {p_{Y_{\kappa}}(y)} \right.\sim\frac{^{{- y}/2}}{\sqrt{2\pi \; y}}} = {{p_{\chi^{2}{(1)}}(y)}.}$

In other words, the PDF of Y_(κ) converges (pointwise) to the PDF of achi-square random variable with one degree of freedom as κ→∞. It canalso be shown that the cumulative distribution function (CDF) of Y_(κ)converges uniformly to the CDF of χ² (1) as κ→∞. Therefore, for large κ,Y_(κ){dot over (˜)}χ²(1) and hence M{dot over (˜)}χ²(n+1).

In principle, one can derive the exact PDF of M using the standardchange of variables theorem in conjunction with the PDFs of thechi-square distribution and (27), though an analytic expression for theresulting convolution is intractable. In practice and in the context ofthe present invention, κ is sufficiently large so that the approximationof Y_(κ) by a χ²(1) random variable (and hence M by χ²(n+1)) isjustified.

III.I.i. Interpretation

The interpretation of the Mahalanobis von Mises statistic (26) can beunderstood by first considering the more general setting of amultivariate random vector x with support on a differentiable manifold

and a PDF of the form p(x)=e^(−ƒ(x)/2). Suppose a point x_(*)ε

is given and one wishes to test the null hypothesis H₀ that x_(*) is nota statistically significant realization of the random vector x (i.e.,x_(*) is a representative draw from x). The p-value for a one-sided testis

P _(Y) _(κ) (y)=p=Pr[ƒ(x)εΩ_(*)]=∫_(Ω*) e ^(−ƒ() x)/2 dx,  (28)

where Ω_(*)={x|ƒ(x)>ƒ(x_(*))≡C}. Smaller p-values imply that therealization x_(*) lies farther out on the tails of the PDF (see FIG. 8a). The null hypothesis H₀ is rejected at the significance level α(typically 0.05 or 0.01) if p<α. FIG. 8 b shows the setup for theanalogous two-sided hypothesis test. For a given significance level α,one determines the contours C_(L) and C_(U) such that

${{\int_{\Omega_{L}}{^{{- {f{(x)}}}/2}{x}}} = {{\int_{\Omega_{U}}{^{{- {f{(x)}}}/2}{x}}} = {\frac{1}{2}\alpha}}},$

where Ω_(L)={x|ƒ(x)<C_(L)} and Ω_(U)={x|ƒ(x)>C_(U)}. (Note that theshaded region in FIG. 8 b has probability α.) A two-sided test withsignificance level a rejects the null hypothesis H₀ if ƒ(x_(*))<C_(L) orƒ(x_(*))>C_(U).

Specializing the one-sided statistical significance test to the Gaussvon Mises distribution, suppose (x,θ)˜GVM(μ, P, α, β, Γ, κ) and arealization (x_(*),θ_(*))ε

^(n)×

is given. For this case,

$\quad\begin{matrix}{{f\left( {x,\theta} \right)} = {{- 2}\ln \; {{\mathcal{M}}\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right)}}} \\{{= {{M\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right)} + {\ln \mspace{14mu} {\det \left( {2\; \pi \; P} \right)}} + {2\; {\ln \left( {2\; \pi \; ^{- \kappa}{I_{0}(\kappa)}} \right)}}}},}\end{matrix}$

where M is the Mahalanobis von Mises statistic ((26)). The integrationregion Ω_(*) in (28) is

$\quad\begin{matrix}{\Omega_{*} = \left\{ {\left( {x,\theta} \right) \in {{\mathbb{R}}^{n} \times }} \middle| {{f\left( {x,\theta} \right)} > {f\left( {x_{*},\theta_{*}} \right)}} \right\}} \\{= {\begin{Bmatrix}\left. {\left( {x,\theta} \right) \in {{\mathbb{R}}^{n} \times }} \middle| {{M\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right)} >} \right. \\{M\left( {x_{*},{\theta_{*};\mu},P,\alpha,\beta,\Gamma,\kappa} \right)}\end{Bmatrix}.}}\end{matrix}$

Substituting this information into (28) yields

$\quad\begin{matrix}{p = {\Pr \left\lbrack {{f\left( {x,\theta} \right)} \in \Omega_{*}} \right\rbrack}} \\{= {\Pr \left\lbrack {{M\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right)} > {M\left( {x_{*},{\theta_{*};\mu},P,\alpha,\beta,\Gamma,\kappa} \right)}} \right\rbrack}} \\{= {\Pr \left\lbrack {{{\chi^{2}(n)} + Y_{\kappa}} > {M\left( {x_{*},{\theta_{*};\mu},P,\alpha,\beta,\Gamma,\kappa} \right)}} \right\rbrack}} \\{\approx {{\Pr \left\lbrack {{\chi^{2}\left( {n + 1} \right)} > {M\left( {x_{*},{\theta_{*};\mu},P,\alpha,\beta,\Gamma,\kappa} \right)}} \right\rbrack}.}}\end{matrix}$

The justification of the approximation in the last line is given in thediscussion proceeding Equation (27). Thus, the resulting p-value can becomputed by evaluating the complementary cumulative distributionfunction (i.e., tail distribution) of χ²(n+1) at the Mahalanobis vonMises statistic (26) evaluated at the realization (x_(*),θ_(*)).

III.J. Asymptotic Expansions

In order to avoid numerical overflow issues, large κ asymptoticexpansions can be used for some expressions containing exponentials andmodified Bessel functions in κ. In what follows, expansions are accurateto 16 significant digits for κ≧500.

For the normalization constant appearing in the definition of the vonMises PDF in (8),

$\begin{matrix}{{\ln \left( {2\pi \; ^{- \kappa}{I_{0}(\kappa)}} \right)} \approx {{\frac{1}{2}{\ln \left( \frac{2\pi}{\kappa} \right)}} + {\frac{1}{8\; \kappa}{\left( {1 + \frac{1}{2\; \kappa} + \frac{25}{48\; \kappa^{2}} + \frac{13}{16\; \kappa^{3}} + \frac{1073}{640\; \kappa^{4}}} \right).}}}} & (29)\end{matrix}$

For the κ-dependent terms in the expression of the differential entropyin (17),

$\begin{matrix}{{{\ln \left( {2\pi \; {I_{0}(\kappa)}} \right)} - {\kappa \; \frac{I_{1}(\kappa)}{I_{0}(\kappa)}}} \approx {{\frac{1}{2}{\ln \left( \frac{2\pi \; }{\kappa} \right)}} + {\frac{1}{4\; \kappa}{\left( {1 + \frac{3}{4\; \kappa} + \frac{25}{24\; \kappa^{2}} + \frac{65}{32\; \kappa^{3}} + \frac{3219}{640\kappa^{4}}} \right).}}}} & (30)\end{matrix}$

For the weights and nodes used in the Gauss von Mises quadratureformulas (derived in Section IV),

$\quad\begin{matrix}{\frac{{B_{1}(\kappa)}^{2}}{{4\; {B_{1}(\kappa)}} - {B_{2}(\kappa)}} \approx \begin{matrix}{{\frac{1}{6} - {\frac{1}{48\kappa^{2}}\begin{pmatrix}{1 + \frac{3}{\kappa} + \frac{151}{16\; \kappa^{2}} +} \\{\frac{267}{8\; \kappa^{3}} + \frac{16967}{128\; \kappa^{4}}}\end{pmatrix}}},{\arccos \left( {\frac{B_{2}(\kappa)}{2\; {B_{1}(\kappa)}} - 1} \right)}} & (31)\end{matrix}} \\{\approx \begin{matrix}{{\sqrt{\frac{3}{\kappa}}\begin{pmatrix}{1 + \frac{1}{4\; \kappa} + \frac{43}{160\; \kappa^{2}} + \frac{443}{896\; \kappa^{3}} + \frac{2523}{2048\; \kappa^{4}} +} \\{\frac{340453}{90112\; \kappa^{5}} + \frac{11584095}{851968\; \kappa^{6}}}\end{pmatrix}},} & (32)\end{matrix}}\end{matrix}$

IV. Gauss von Mises Quadrature

The Julier-Uhlmann unscented transform or the more general framework ofGauss-Hermite quadrature enables one to compute the expected value of anon-linear transformation of a multivariate Gaussian random vector. Inthis section, these methodologies are extended to enable the computationof the expected value a non-linear transformation of a Gauss von Mises(GVM) random vector, thereby providing a general framework for Gauss vonMises quadrature. The quadrature formulas are subsequently used in theprediction step of the GVM filter developed in the next section. The GVMquadrature framework is also useful for extracting supplementalstatistics from a GVM distribution. For example, if a space object'sorbital state is represented as a GVM distribution in equinoctialorbital element space, one can compute the expected value of theobject's Cartesian position and velocity or the covariance in itsposition and velocity.

Given (x,θ)˜GVM(μ, P, α, β, Γ, κ) and a function ƒ:

^(n)×

→

, an approximation is sought for

E[ƒ(x,θ)]=

_(n)∫_(−π) ^(π)

(x,θ;μ,P,α,β,Γ,κ)ƒ(x,θ)dθdx.  (33)

As in classical Gaussian quadrature, the framework of Gauss von Misesquadrature approximates (33) as a weighted sum of function values atspecified points:

$\begin{matrix}{{\int_{{\mathbb{R}}^{n}}{\int_{- \pi}^{\pi}{\; \; {\mathcal{M}\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right)}{f\left( {x,\theta} \right)}{\theta}{x}}}} \approx {\sum\limits_{i = 1}^{N}{w_{\sigma_{i}}{{f\left( {x_{\sigma_{i}},\theta_{\sigma_{i}}} \right)}.}}}} & (34)\end{matrix}$

The set of quadrature nodes, sometimes called sigma points, {(x_(σ) _(i),θ_(σ) _(i) )}_(i=1) ^(N) and corresponding quadrature weights {w_(σ)_(i) }_(i=1) ^(N) can be chosen so that (34) is exact for a certainclass of functions. In the derivation of the GVM quadrature weights andnodes, the Smolyak sparse grid paradigm can be used so that, for aspecified order of accuracy, the number of nodes increases polynomiallywith the dimension n thereby avoiding the so-called “curse ofdimensionality.” In particular, the number of quadrature nodes (andhence function evaluations) increases linearly in the dimension n forthe third-order rule derived in Subsection IV.A and increasesquadratically in n for the fifth-order rule derived in Subsection IV.B.Higher-order GVM quadrature rules are discussed in Subsection IV.C.

IV.A. Third-Order Method

Without loss of generality, the method can be restricted to quadratureusing the canonical GVM distribution as the weighting function:

$\begin{matrix}{{\int_{{\mathbb{R}}^{n}}{\int_{- \pi}^{\pi}{\left( {{z;0},1} \right)\; {\mathcal{M}\left( {{\varphi;0},\kappa} \right)}{f\left( {z,\varphi} \right)}{\varphi}{z}}}} \approx {\sum\limits_{i = 1}^{N}{w_{\sigma_{i}}{{f\left( {z_{\sigma_{i}},\varphi_{\sigma_{i}}} \right)}.}}}} & (35)\end{matrix}$

Indeed, if considering the general quadrature problem in Equation (34),the nodes (x_(σ) _(i) , θ_(σ) _(i) ) can be recovered from the nodes(z_(σ) _(i) , φ_(σ) _(i) ) of the canonical problem (35) through thetransformation (20):

${x_{\sigma_{i}} = {\mu + {A\; z_{\sigma_{i}}}}},{\theta_{\sigma_{i}} = {\varphi_{\sigma_{i}} + \alpha + {\beta^{T}z_{\sigma_{i}}} + {\frac{1}{2}z_{\sigma_{i}}^{T}\Gamma \; {z_{\sigma_{i}}.}}}}$

In what follows, the following notation is used. Let ξε

, ηε(−π,π], and

={(i,j)ε

²|1≦i<j≦n},

⁰⁰={(z,φ)ε

^(n)×

|z=0,φ=0},

^(η0)={(z,φ)ε

^(n)×

|z=0,|φ|=η},

_(i) ^(ξ0)={(z,φ)ε

^(n)×

|z₁ = . . . =z _(i−1) =z _(i+1) = . . . =z _(n)=0,|z _(i)|=ξ,φ=0},

_(i) ^(ξη)={(z,φ)ε

^(n)×

|z₁ = . . . =z _(i−1) =z _(i+1) = . . . =z _(n)=0,|z _(i)|=ξ,|φ|=η},

_(ij) ^(ξξ)={(z,φ)ε

^(n)×

|z₁ = . . . =z _(i−1) =z _(i+1) = . . . =z _(j−1) =z _(j+1) = . . . =z_(n)=0,

|z _(i) |=|z _(j)|=ξ,φ=0}.

For the third-order method, it is demanded that the approximation (35)be exact for all functions of the form

$\begin{matrix}{{{f\left( {z,\varphi} \right)} = {\left( {a_{0} + {\frac{1}{2}z^{T}B_{0}z}} \right) + {a_{1}\cos \; \varphi} + {a_{2}\cos \; 2\varphi} + {g_{o}\left( {z,\varphi} \right)}}},} & (36)\end{matrix}$

where a₀, a₁, a₂ε

, B₀ is any symmetric n×n matrix, and g_(o) is any function such thatg_(o)(z,φ)=−g_(o)(−z,φ) or g_(o)(z,φ)=−g(z,−φ). It is claimed that thisobjective can be achieved with the quadrature node set

= 00 ⋃ η   0 ⋃ ⋃ n i = 1  i ξ   0 , ( 37 )

for some choice of parameters ξ and η, and a quadrature rule of the form

∫ ℝ n  ∫ - π π    ( z ; 0 , I )   ( φ ; 0 , κ )  f  ( z , φ )  φ    z ≈ w 00  ∑ ( z , φ ) ∈ 00   f  ( z , φ ) + w η0  ∑ ( z ,φ ) ∈ η0   f  ( z , φ ) + w ξ0  ∑ i = 1 n   ∑ ( z , φ ) ∈ i ξ0  f  ( z , φ ) . ( 38 )

Indeed, to make (38) exact for functions of the form (36), it sufficesto impose the condition that (38) be exact for the basis functionslisted in Table 1.

TABLE 1 f (z, φ) Constraint 1 w₀₀ + 2w_(η0) + 2nw_(ξ0) = 1 cos φ w₀₀ +2cos η w_(η0) + 2nw_(ξ0) = I₁(K)/I₀(K) cos 2φ w₀₀ + 2cos2η w_(η0) +2nw_(ξ0) = I₂(K)/I₀(K) z_(i) ², i ∈ {1, . . . , n} 2ξ²w_(ξ0) = 1 z_(i)⁴, i ∈ {1, . . . , n} 2ξ⁴w_(ξ0) = 3The corresponding constraints on the quadrature weights and theparameters ξ and η appearing in the quadrature nodes are also listed inthe table. It is noted that these constraints follow from the assumedquadrature rule (38) in conjunction with the characteristic function(16) or the low-order moments listed in Section III.C. Solving the fiveconstraint equations in Table 1 for w₀₀, w_(η0), w_(ξ0), ξ, and η yields

$\begin{matrix}{{\xi = \sqrt{3}},{\eta = {\arccos \left( {\frac{B_{2}(\kappa)}{2{B_{1}(\kappa)}} - 1} \right)}},{w_{\xi 0} = \frac{1}{6}},{w_{\eta 0} = \frac{{B_{1}(\kappa)}^{2}}{{4{B_{1}(\kappa)}} - {B_{2}(\kappa)}}},{w_{00} = {1 - {2w_{\eta 0}} - {2{nw}_{\xi 0}}}},} & (39)\end{matrix}$

where B_(p)(κ)≡1−I_(p)(κ)/I₀(κ). For large κ, the asymptotic expansionsin (31) and (32) should be applied to compute the quantities containingthe modified Bessel functions.

In summary, the third-order GVM quadrature rule is (38) with the nodesand weights specified in (37) and (39).

It is noted that the number of quadrature nodes in the set (37) is

|

|=2n+3.

Thus, the number of nodes increases linearly with the dimension n.Further, this is precisely the same number of nodes (sigma points) usedin the unscented transform (which is a third-order Gauss-Hermitequadrature method) in a Cartesian space of dimension n+1.

IV.B. Fifth-Order Method

The derivation of the fifth-order GVM quadrature rule is similar to thatof the third-order rule. It is demanded that the approximation (20) beexact for all functions of the form

$\begin{matrix}{{{f\left( {z,\varphi} \right)} = {\left( {a_{0} + {\frac{1}{2}b_{0}^{ij}z_{i}z_{j}} + {\frac{1}{24}c_{0}^{{ijk}\; }z_{i}z_{i}z_{k}z_{}}} \right) + {\left( {a_{1} + {\frac{1}{2}b_{1}^{ij}z_{i}z_{j}}} \right)\cos \; \varphi} + {a_{2}\cos \; 2\varphi} + {g_{o}\left( {z,\varphi} \right)}}},} & (40)\end{matrix}$

where the coefficients are arbitrary and g_(o) is any function such thatg_(o)(z,φ)=−g_(o)(−z,φ) or g_(o)(z,φ)=−g(z,−φ). (The summationconvention is implied between repeated upper and lower indices.) It isclaimed that this objective can be achieved with the quadrature node set

$\begin{matrix}{{ = {^{00}\bigcup ^{\eta 0}\bigcup{\bigcup\limits_{i = 1}^{n}_{i}^{\xi 0}}\bigcup{\bigcup\limits_{i = 1}^{n}_{i}^{\xi\eta}}\bigcup{\bigcup\limits_{{({i,j})} \in }_{ij}^{\xi\xi}}}},} & (41)\end{matrix}$

for some choice of parameters ξ and η, and a quadrature rule of the form

∫ ℝ n  ∫ - π π    ( z ; 0 , I )   ( φ ; 0 , κ )  f  ( z , φ )  φ    z ≈ w 00  ∑ ( z , φ ) ∈ 00   f  ( z , φ ) + w η0  ∑ ( z ,φ ) ∈ η0   f  ( z , φ ) + w ξ0  ∑ i = 1 n  ∑ ( z , φ ) ∈ i ξ0   f ( z , φ ) + w ξη  ∑ i = 1 n  ∑ ( z , φ ) ∈ i ξη  f  ( z , φ ) + wξξ  ∑ ( i , j ) ∈    ∑ ( z , φ ) ∈ ij ξξ  f  ( z , φ ) .  ( 42 )

Indeed, to make (42) exact for functions of the form (40), it sufficesto impose the condition that (42) be exact for the basis functionslisted in Table 2.

TABLE 2 f (z, φ) Constraint 1${w_{00} + {2w_{\eta \; 0}} + {2{nw}_{\xi \; 0}} + {4{nw}_{\xi\eta}} + {4\begin{pmatrix}n \\2\end{pmatrix}w_{\xi\xi}}} = 1$ cos φ${w_{00} + {2\cos \; \eta \; w_{\eta \; 0}} + {2{nw}_{\xi \; 0}} + {4n\; \cos \; \eta \; w_{\xi\eta}} + {4\begin{pmatrix}n \\2\end{pmatrix}w_{\xi\xi}}} = {{I_{1}(\kappa)}/{I_{0}(\kappa)}}$ cos 2φ${w_{00} + {2\cos \; \eta \; w_{\eta \; 0}} + {2{nw}_{\xi \; 0}} + {4n\; \cos \; \eta \; w_{\xi\eta}} + {4\begin{pmatrix}n \\2\end{pmatrix}w_{\xi\xi}}} = {{I_{2}(\kappa)}/{I_{0}(\kappa)}}$ z_(i)², i ∈ {1, . . . , n} 2ξ²w_(ξ0) + 4ξ²w_(ξη) + 4(n − 1)ξ²w_(ξξ) = 1 z_(i)⁴, i ∈ {1, . . . , n} 2ξ⁴w_(ξ0) + 4ξ⁴w_(ξη) + 4(n − 1)ξ⁴w_(ξξ) = 3 z_(i)²z_(j) ², (i, j) ∈ P 4ξ⁴w_(ξξ) = 1 z_(i) ²cos φ, i ∈ {1, . . . , n}2ξ²w_(ξ0) + 4ξ²cosηw_(ξη) + 4(n − 1)ξ²w_(ξξ) = I₁(κ)/I₀(κ)The corresponding constraints on the quadrature weights and theparameters ξ and η appearing in the quadrature nodes are also listed inthe table. Solving the seven constraint equations in Table 2 for w₀₀,w_(η0), w_(ξ0), w_(ξη), w_(ξξ), ξ, and η yields

$\begin{matrix}{{{\xi = \sqrt{3}},{\eta = {\arccos \left( {\frac{B_{2}(\kappa)}{2{B_{1}(\kappa)}} - 1} \right)}},{w_{\xi\xi} = \frac{1}{36}},{w_{\xi\eta} = {\frac{1}{6} \cdot \frac{{B_{1}(\kappa)}^{2}}{{4{B_{1}(\kappa)}} - {B_{2}(\kappa)}}}},{w_{\xi 0} = {\frac{2}{9} - \frac{n}{18} - {\frac{1}{3} \cdot \frac{{B_{1}(\kappa)}^{2}}{{4{B_{1}(\kappa)}} - {B_{2}(\kappa)}}}}},{w_{\eta 0} = {\left( {1 - \frac{n}{3}} \right) \cdot \frac{{B_{1}(\kappa)}^{2}}{{4{B_{1}(\kappa)}} - {B_{2}(\kappa)}}}},{w_{00} = {1 - {2w_{\eta 0}} - {2{nw}_{\xi 0}} - {4{nw}_{\xi\eta}} - {2{n\left( {n - 1} \right)}w_{\xi\xi}}}},{where}}{{B_{p}(\kappa)} \equiv {1 - {{I_{p}(\kappa)}/{{I_{0}(\kappa)}.}}}}} & (43)\end{matrix}$

In summary, the fifth-order GVM quadrature rule is (42) with the nodesand weights specified in (41) and (43).

It is noted that the number of quadrature nodes in the set (41) is

$\begin{matrix}{{} = {1 + 2 + {2n} + {4n} + {4\begin{pmatrix}n \\2\end{pmatrix}}}} \\{= {{2\left( {n + 1} \right)^{2}} + 1.}}\end{matrix}$

Thus, the number of nodes increases quadratically with the dimension n.Further, this is precisely the same number of nodes obtained using thesparse grid fifth-order Gauss-Hermite quadrature method of Genz andKeister in a Cartesian space of dimension n+1.

IV.C. Higher-Order Methods

It is briefly noted that that GVM quadrature rules can be generated toarbitrarily high-order using classical Gauss-Hermite quadrature formulasand the discrete Fourier transform. Indeed,

${{ \equiv {\int_{{\mathbb{R}}^{n}}^{\;}{\int_{- \pi}^{\pi}{\ \left( {{z;0},I} \right)\left( {{\varphi;0},\kappa} \right){f\left( {z,\varphi} \right)}{\varphi}\ {z}}}}} = {{\int_{- \pi}^{\pi}{\left( {{\varphi;0},\kappa} \right){\int_{{\mathbb{R}}^{n}}^{\;}{\ \left( {{z;0},I} \right)f\left( {z,\varphi} \right){z}\ {\varphi}}}}} \approx {\sum\limits_{j = 1}^{N}\; {w_{\sigma \; j}^{GH}{\int_{- \pi}^{\pi}{{\left( {{\varphi;0},\kappa} \right)}{f\left( {z_{\sigma_{j}}^{GH},\varphi} \right)}{\varphi}}}}}}},$

where {w_(σj) ^(GH)}_(j=1) ^(N) and {z_(σj) ^(GH)}_(j=1) ^(N) are thecanonical n-dimensional Gauss-Hermite weights and nodes. Thesequadrature weights and nodes can be generated using, for example, thesparse grid framework of Genz and Keister, to yield an approximation toany desired order of accuracy. For each j=0, . . . , N, perform adiscrete Fourier transform on ƒ(z_(σj) ^(GH),φ) and introduce thenotation

$\mspace{76mu} {{{f\left( {z_{\sigma_{j}}^{GH},\varphi} \right)} \approx {\sum\limits_{k = {- M}}^{M}\; {F_{jk}^{\; k\; \varphi}}}},\mspace{79mu} {where}}$$\mspace{79mu} {{F_{jk} = {\frac{1}{{2M} + 1}{\sum\limits_{ = {- M}}^{M}\; {{f\left( {z_{\sigma_{j}}^{GH},\varphi_{}} \right)}^{- {k\varphi}_{}}}}}},\mspace{79mu} {\varphi_{} = {\frac{2{\pi }}{{2M} + 1}.\mspace{79mu} {Therefore}}},{{ \approx {\sum\limits_{j = 1}^{N}\; {w_{\sigma_{j}}^{GH}{\sum\limits_{k = {- M}}^{M}\; {F_{jk}{\int_{- \pi}^{\pi}{\left( {{\varphi;0},\kappa} \right)^{\; k\; \varphi}\ {\varphi}}}}}}}} = {\sum\limits_{j = 1}^{N}\; {w_{\sigma_{j}}^{GH}{\sum\limits_{k = {- M}}^{M}\; {\frac{I_{k}(\kappa)}{I_{0}(\kappa)}{F_{jk}.}}}}}}}$

In the context of the present invention which makes use of GVMquadrature, it is observed that there is little benefit in using anyquadrature rule beyond the simple third-order method.

V. Uncertainty Propagation

This section presents the method to implement the prediction step of theBayesian non-linear filter using the Gauss von Mises (GVM) probabilitydensity function (PDF) as input. Effectively, the new algorithm providesa means for approximating the non-linear transformation of a GVMdistribution as another GVM distribution. The GVM quadrature rulesdeveloped in the previous section play a key role in the computation. Inanalogy to the unscented transform applicable to Gaussian PDFs,quadrature nodes can be deterministically selected from the initial GVMdistribution and then acted on by the non-linear transformation. Thetransformed quadrature nodes can then be used to reconstruct theparameters of the transformed GVM distribution. One application of thismethodology can be the propagation of a space object's state uncertaintyunder non-linear two-body dynamics, which provides improved predictioncapabilities of the object's future location and characterization of itsorbital uncertainty at future times. It is shown in the EXAMPLE thatuncertainty propagation using the new GVM filter prediction step forthis application can be achieved by propagating 13 quadrature nodes, thesame number used in the standard unscented Kalman filter (UKF).Moreover, the new method is shown to maintain a proper characterizationof the orbital state uncertainty for up to eight times as long as theUKF.

The organization of this section is as follows. In Subsection V.A, thepreliminary notation is defined followed by discussions on how statepropagation under a non-linear system of ordinary differential equations(ODEs) fits into the general framework. In Subsection V.B, the GVMfilter prediction step is motivated followed by the complete algorithmdescription in Subsection V.C. In Subsection V.D, two different onlinemetrics are proposed which allow the operator to validate thecharacterization of the transformed PDF by a GVM distribution. InSubsection V.E, it is shown how the inclusion of additional uncertainparameters or stochastic process noise can be treated within the sameframework. Finally, in Subsection V.F, it is shown how a “mixtureversion” of the GVM filter prediction step can be formulated, in analogyto the Gaussian sum (mixture) filter to provide proper uncertaintyrealism in the most challenging scenarios.

V.A. Preliminary Notation

Let (x,θ)˜GVM(μ, P, α, β, Γ, κ), Φ:

^(n)×

→

^(n)×

be a diffeomorphism (i.e., a bijection on

^(n)×

with a smooth inverse), and ({tilde over (x)}, {tilde over (θ)}) bedefined such that

({tilde over (x)},{tilde over (θ)})=Φ(x,θ;p),  (44)

where p denotes any constant non-stochastic parameters. The inverse of(44) is denoted as

(x,θ)=Ψ({tilde over (x)},{tilde over (θ)};p).  (45)

Thus, given the GVM random vector (x,θ) and a diffeomorphism Φ, theobjective is to approximate the joint PDF of ({tilde over (x)},{tildeover (θ)}) by a GVM distribution and quantify the fidelity of thisapproximation.

The first-order system of ODEs,

′(t)=

(

(t),t),  (46)

naturally gives rise to a family of diffeomorphisms. For a given initialcondition

(t₀)=x₀, the solution of (46) is denoted as

(t)=Φ(

₀ ;t ₀ ,t),  (47)

which maps the state

₀ at some initial epoch to the state at a future time. (Existence anduniqueness of solutions is assumed on the interval [t₀,t].) In thiscontext, Φ is called the solution flow and is of the form (44) where theparameter vector p contains the initial and final times. The inversesolution flow is denoted as

₀=Ψ(x(t);t ₀ ,t),

which maps the state

(t) at time t to the state x₀ at some past time t₀. If the initialcondition x₀ is uncertain and described by a PDF, then Equation (47)dictates how this PDF is transformed to a future epoch.

The dynamics governing the two-body problem in orbital mechanics are ofthe form (46) where

encodes the force components (e.g., gravity, atmospheric drag, etc.)and, typically, the state

is six-dimensional Cartesian Earth Centered Inertial (ECI)position-velocity coordinates. (The choice of inertial coordinates canbe made so that Newton's laws hold.) As argued in Section I.A, it isadvantageous to represent orbital states and uncertainties inequinoctial orbital element (EqOE) coordinates because the solution flow(47) in such coordinates, denoted as Φ^(EqOE), is approximated moreclosely as a transformation of the form (18) which maps a GVMdistribution to a GVM distribution. Computationally, one can implementthe map Φ^(EqOE) by transforming the ODEs expressed in the natural ECIcoordinates to a system of ODEs in EqOE space (see Equation (2)) andthus propagating states directly in EqOE space. In an equivalentapproach, one can (i) take the initial condition

₀ ^(EqOE) in EqOE space and convert it to an initial condition

₀ ^(ECI) in ECI coordinates, (ii) propagate this initial condition usingthe solution flow Φ^(ECI) with respect to ECI coordinates to obtain thepropagated ECI state

_(ƒ) ^(ECI), and (iii) convert

_(ƒ) ^(ECI) back to EqOE coordinates to yield the propagated EqOE state

_(ƒ) ^(EqOE). Thus, Φ^(EqOE) is a composition of the three mapsdescribed in Steps (i)-(iii). This approach is sometimes preferred whenusing a commercial orbital propagator utilizing input in ECIcoordinates. In summary, the following diagram commutes:

In what follows, the following notation is used. Given (x,θ)ε

^(n)×

and a diffeomorphism Φ:

^(n)×

→

^(n)×

, the following definitions are used.

-   -   Φ_(x):        ^(n)×        →        ^(n) is the projection of Φ on        ^(n).    -   Φ_(θ):        ^(n)×        →        is the projection of Φ on        .    -   ∂_(x)Φ_(x):        ^(n)×        →        ^(n×n) is the Jacobian of Φ_(x) with respect to x:

${\partial_{x}{\Phi_{x}\left( {x,\theta} \right)}} = \begin{bmatrix}\frac{\partial{\Phi_{1}\left( {\xi,\theta} \right)}}{\partial\xi_{1}} & \ldots & \frac{\partial{\Phi_{1}\left( {\xi,\theta} \right)}}{\partial\xi_{n}} \\\vdots & \ddots & \vdots \\\frac{\partial{\Phi_{n}\left( {\xi,\theta} \right)}}{\partial\xi_{1}} & \ldots & \frac{\partial{\Phi_{n}\left( {\xi,\theta} \right)}}{\partial\xi_{n}}\end{bmatrix}_{\xi = x}$

-   -   ∂_(x)Φ_(θ):        ^(n)×        →        ^(n) is the gradient of Φ_(θ) with respect to x:

${\partial_{x}{\Phi_{\theta}\left( {x,\theta} \right)}} = \begin{bmatrix}\frac{\partial{\Phi_{\theta}\left( {\xi,\theta} \right)}}{\partial\xi_{1}} \\\vdots \\\frac{\partial{\Phi_{\theta}\left( {\xi,\theta} \right)}}{\partial\xi_{n}}\end{bmatrix}_{\xi = x}$

-   -   ∂_(x) ²Φ_(θ):        ^(n)×        →        ^(n×n) is the Hessian of Φ_(θ) with respect to x:

${\partial_{x}^{2}{\Phi_{\theta}\left( {x,\theta} \right)}} = \begin{bmatrix}\frac{\partial^{2}{\Phi_{\theta}\left( {\xi,\theta} \right)}}{{\partial\xi_{1}}{\partial\xi_{1}}} & \ldots & \frac{\partial^{2}{\Phi_{\theta}\left( {\xi,\theta} \right)}}{{\partial\xi_{1}}{\partial\xi_{n}}} \\\vdots & \ddots & \vdots \\\frac{\partial^{2}{\Phi_{\theta}\left( {\xi,\theta} \right)}}{{\partial\xi_{n}}{\partial\xi_{1}}} & \ldots & \frac{\partial^{2}{\Phi_{\theta}\left( {\xi,\theta} \right)}}{{\partial\xi_{n}}{\partial\xi_{n}}}\end{bmatrix}_{\xi = x}$

V.B. Motivation

Consider the random variables ({tilde over (x)},{tilde over (θ)})defined by

({tilde over (x)},{tilde over (θ)})=Φ(x,θ;p),

where Φ:

^(n)×

→

^(n)×

is a diffeomorphism and (x,θ)˜GVM(μ, P, α, β, Γ, κ). By the change ofvariables theorem for PDFs, it follows that

$\begin{matrix}{{{p\left( {\overset{\sim}{x},\overset{\sim}{\theta}} \right)} = {{\det \left\lbrack \frac{\partial\Psi}{\partial x} \right\rbrack}\left( {{\Psi_{x}\left( {\overset{\sim}{x},{\overset{\sim}{\theta};p}} \right)},{{\Psi_{\theta}\left( {\overset{\sim}{x},{\overset{\sim}{\theta};p}} \right)};\mu},P,\alpha,\beta,\Gamma,\kappa} \right)}},} & (49)\end{matrix}$

where Ψ is the inverse of Φ and {tilde over (

)}=({tilde over (x)},{tilde over (θ)}). The transformed PDF (49) is aGVM distribution in ({tilde over (x)},{tilde over (θ)}) if thetransformation Φ (or Ψ) is of the form (18). In general, the non-lineartransformation of a GVM distribution is not a GVM distribution. In suchcases, the objective of the GVM filter prediction step is to compute anapproximate GVM distribution such that

p({tilde over (x)},{tilde over (θ)})≈

({tilde over (x)},{tilde over (θ)};{tilde over (μ)}{tilde over(P)},{tilde over (α)},{tilde over (β)},{tilde over (Γ)},{tilde over(κ)}).  (50)

One way to achieve this objective is to approximate Φ as a Taylor seriesabout the modal point (x,θ)=(μ,α) and then match the Taylor coefficientswith the coefficients in the transformation (18). It follows that

$\begin{matrix}{{L = {\partial_{x}{\Phi_{x}\left( {\mu,\alpha} \right)}}},{d = {{\Phi_{x}\left( {\mu,\alpha} \right)} - {{\partial_{x}{\Phi_{x}\left( {\mu,\alpha} \right)}}\mu}}},{a = {{\Phi_{\theta}\left( {\mu,\alpha} \right)} - \alpha - {{\partial_{x}{\Phi_{\theta}\left( {\mu,\alpha} \right)}^{T}}\mu} + {\frac{1}{2}\mu^{T}{\partial_{x}^{2}{\Phi_{\theta}\left( {\mu,\alpha} \right)}}\mu}}},{b = {{\partial_{x}{\Phi_{\theta}\left( {\mu,\alpha} \right)}} - {{\partial_{x}^{2}{\Phi_{\theta}\left( {\mu,\alpha} \right)}}\mu}}},{C = {\partial_{x}^{2}{\Phi_{\theta}\left( {\mu,\alpha} \right)}}},} & (51)\end{matrix}$

where, for notational convenience, the dependence of Φ and itsderivatives on the non-stochastic parameter p is suppressed.Substituting (51) into (19) yields

{tilde over (μ)}=Φ_(x)(μ,α),{tilde over (P)}=∂ _(x)Φ_(x)(μ,α)P∂_(x)Φ_(x)(μ,α)^(T),{tilde over (α)}=Φ_(θ)(μ,α),

{tilde over (β)}=Q^(T) [β+A ^(T)∂_(x)Φ_(θ)(μ,α)],{tilde over (Γ)}=Q^(T)[Γ+A ^(T)∂_(x) ²Φ_(θ)(μ,α)A]Q,{tilde over (κ)}=κ,  (52)

where Q is an n×n orthogonal matrix and R is an n×n positive-definiteupper-triangular matrix such that QR=(LA)^(T). It follows thatQ=A^(T)∂_(x)Φ_(x)(μ,α)^(T)Ã^(−T) where Ã is the lower-triangular matrixsuch that {tilde over (P)}=ÃÃ^(T); hence

{tilde over (μ)}=Φ_(x)(μ,α),{tilde over (P)}=∂ _(x)Φ_(x)(μ,α)P∂_(x)Φ_(x)(μ,α)^(T),{tilde over (α)}=Φ_(θ)(μ,α),

{tilde over (β)}=Ã⁻¹∂_(x)Φ_(x)(μ,α)A[β+A ^(T)∂_(x)Φ_(θ)(μ,α)],

{tilde over (Γ)}=Ã⁻¹∂_(x)Φ_(x)(μ,α)A[Γ+A ^(T)∂_(x) ²Φ_(θ)(μ,α)A]A^(T)∂_(x)Φ_(x)(μ,α)^(T) Ã ^(−T),{tilde over (κ)}=κ.  (53)

An equivalent form of the above equations for {tilde over (β)} and{tilde over (Γ)} are used in the next section which expresses allpartial derivatives of Φ with respect to the standardized variablez=A⁻¹(x−μ). By the chain rule,

∂_(x)Φ_(θ) =A ^(−T)∂_(z)Φ_(θ),∂_(x) ²Φ_(θ) =A ^(−T)∂_(z) ²Φ_(θ) A⁻¹,∂_(x)Φ_(x)=∂_(z)Φ_(x) A ⁻¹.

Substituting the above equations into (53) yields

{tilde over (μ)}=Φ_(x)(μ,α),{tilde over (P)}=∂ _(x)Φ_(x)(μ,α)P∂_(x)Φ_(x)(μ,α)^(T),{tilde over (α)}=Φ_(θ)(μ,α),

{tilde over (β)}=Ã⁻¹∂_(z)Φ_(x)(μ,α)[β+∂_(z)Φ_(θ)(μ,α)],

{tilde over (Γ)}=Ã⁻¹∂_(z)Φ_(x)(μ,α)[Γ+∂_(z)²Φ_(θ)(μ,α)]∂_(z)Φ_(x)(μ,α)^(T) Ã ^(−T),{tilde over (κ)}=κ.  (54)

It is noted that Equations (53) or (54) are analogous to those used inthe prediction step of the extended Kalman filter (EKF) for thelinearized propagation of a multivariate Gaussian distribution. Indeed,if x˜N(μ,P), Φ:

^(n)→

^(n) then, under weak non-linear assumptions in Φ, {tilde over(x)}=Φ(x){dot over (˜)}N({tilde over (μ)},{tilde over (P)}), where

{tilde over (μ)}=Φ(μ),{tilde over (P)}=∂ _(x)Φ(μ)P∂ _(x)Φ(μ)^(T).

While the analogous “EKF-like” prediction derived in Equations (53)provides an approximate GVM distribution to a non-linear transformationof a GVM random vector, the proposed algorithm developed in the nextsubsection provides a more accurate approximation while avoiding thedirect computation of the partial derivatives of the transformation Φ.

V.C. Algorithm Description

FIG. 1 is a flowchart illustrating an exemplary process for transforminga Gauss von Mises (GVM) distribution under a non-linear transformationand approximating the output as a GVM distribution. In this example,there are two inputs 100: (i) a random vector (x,θ)ε

^(n)×

distributed as a GVM distribution with parameter set (μ, P, α, β, Γ, κ),and (ii) a diffeomorphism Φ:

^(n)×

→

^(n)×

. Where, as described in detail above: (μ,α) is the mode of the GVMdistribution (i.e., the state with the maximum likelihood); P is thecovariance of x; κ quantifies concentration in the angular variable θ; βbeta characterizes correlation between x and θ; and Γ modelshigher-order cumulants (which give the level sets of the GVMdistribution their distinctive “banana” or “boomerang” shape). Theoutput 105 can be an approximation of the distribution of ({tilde over(x)},{tilde over (θ)})=Φ(x,θ) as a GVM distribution with parameter set({tilde over (μ)}, {tilde over (P)}, {tilde over (α)}, {tilde over (β)},{tilde over (Γ)}, {tilde over (κ)}). The algorithm has four main stepssummarized below:

-   -   Generate 101 N quadrature nodes {(x_(σ) _(i) ,θ_(σ) _(i)        )}_(i=1) ^(N) from the input GVM distribution and compute the        transformed quadrature nodes ({tilde over (x)}_(σ) _(i) ,{tilde        over (θ)}_(σ) _(i) )=Φ(x_(σ) _(i) ,θ_(σ) _(i) ), for i=1, . . .        , N;    -   Recover 102 the parameters {tilde over (μ)} and {tilde over (P)}        of the transformed GVM distribution from the transformed        quadrature nodes computed in Step 101 using the appropriate        Gauss von Mises quadrature rule;    -   Compute 103 approximations of {tilde over (α)}, {tilde over        (β)}, and {tilde over (Γ)}, denoted as {circumflex over (α)},        {circumflex over (β)}, and {circumflex over (Γ)}, using        Equations (54) in conjunction with suitable approximations of        the partial derivatives of Φ; and    -   Set 104 {tilde over (κ)}=κ and recover the parameters {tilde        over (α)}, {tilde over (β)}, and {tilde over (Γ)} by solving the        non-linear least squares problem

$\begin{matrix}{\mspace{734mu} (55)} & \; \\{{\left( {\overset{\sim}{\alpha},\overset{\sim}{\beta},\overset{\sim}{\Gamma}} \right) = {\underset{\hat{\alpha},\hat{\beta},\hat{\Gamma}}{argmin}{\sum\limits_{i = 1}^{N}\left\lbrack {{M\left( {x_{\sigma_{i}},{\theta_{\sigma_{i}};\mu},P,\alpha,\beta,\Gamma,\kappa} \right)} - {M\left( {{\overset{\sim}{x}}_{\sigma_{i}},{{\overset{\sim}{\theta}}_{\sigma_{i}};\overset{\sim}{\mu}},\overset{\sim}{P},\hat{\alpha},\hat{\beta},\hat{\Gamma},\overset{\sim}{\kappa}} \right)}} \right\rbrack^{2}}}},} & \;\end{matrix}$

-   -   -   where M is the Mahalanobis von Mises statistic (26).

These four steps 101, 102, 103, and 104 are described in more detailbelow including discussions on how they can be specialized to the casewhen Φ is the solution flow (in equinoctial orbital element coordinates)corresponding to the two-body problem in orbital mechanics (henceforthreferred to as the “two-body problem”).

At step 101, using either the third, fifth, or higher-order GVMquadrature rule derived in Sections IV.A, IV.B, and IV.C, respectively,first generate N quadrature nodes {(z_(σ) _(i) ,φ_(σ) _(i) )}_(i=1) ^(N)with respective weights {w_(σ) _(i) }_(i=1) ^(N) corresponding to thecanonical GVM distribution. Next, use these canonical quadrature nodesto generate quadrature nodes corresponding to the input GVM randomvector (x,θ)˜GVM(μ, P, α, β, Γ, κ) according to

${x_{\sigma_{i}} = {\mu + {Az}_{\sigma_{i}}}},{\theta_{\sigma_{i}} = {\varphi_{\sigma_{i}} + \alpha + {\beta^{T}z_{\sigma_{i}}} + {\frac{1}{2}z_{\sigma_{i}}^{T}\Gamma \; {z_{\sigma_{i}}.}}}}$

for i=1, . . . , N. Finally, compute the transformed quadrature nodes({tilde over (x)}_(σ) _(i) ,{tilde over (θ)}_(σ) _(i) )=Φ(x_(σ) _(i),θ_(σ) _(i) ), for i=1, . . . , N.

In the case of the two-body problem, the transformed quadrature nodescan be generated by propagating each quadrature node (x_(σ) _(i) ,θ_(σ)_(i) ), representing a state in equinoctial orbital element space, froma specified initial time to a specified final time under the underlyingdynamics. This generally can be accomplished using the numericalsolution of a non-linear system of ordinary differential equations.Additional details are discussed in the paragraph preceding Equation(48).

At step 102, the transformed parameters {tilde over (μ)} and {tilde over(P)} are given by

{tilde over (μ)}=E[{tilde over (x)}]=E[Φ_(x)(x,θ)]

{tilde over (P)}=E[({tilde over (x)}−{tilde over (μ)})({tilde over(x)}−{tilde over (μ)})^(T) ]=E[(Φ_(x)(x,θ)−{tilde over(μ)})(Φ_(x)(x,θ)−{tilde over (μ)})^(T)].

Using the quadrature weights {w_(σ) _(i) }_(i=1) ^(N) and thetransformed quadrature nodes {({tilde over (x)}_(σ) _(i) ,{tilde over(θ)}_(σ) _(i) )}_(i=1) ^(N) generated in Step 101, the above expectedvalues can be approximated using the framework of Gauss von Misesquadrature developed in Section IV:

${\overset{\sim}{\mu} = {\sum\limits_{i = 1}^{N}{w_{\sigma_{i}}{\overset{\sim}{x}}_{\sigma_{i}}}}},{\overset{\sim}{P} = {\sum\limits_{i = 1}^{N}{{w_{\sigma_{i}}\left( {{\overset{\sim}{x}}_{\sigma_{i}} - \overset{\sim}{\mu}} \right)}{\left( {{\overset{\sim}{x}}_{\sigma_{i}} - \overset{\sim}{\mu}} \right)^{T}.}}}}$

The lower-triangular Cholesky factor of {tilde over (P)}, satisfying{tilde over (P)}=ÃÃ^(T), can be used in subsequent computations.

At step 103, the parameters {tilde over (α)}, {tilde over (β)}, and{tilde over (Γ)} are approximated as {circumflex over (α)}, {circumflexover (β)}, and {circumflex over (Γ)} using Equations (54):

{circumflex over (α)}=Φ_(θ)(μ,α),{circumflex over(β)}=Ã⁻¹∂_(z)Φ_(x)(μ,α)[β+∂_(z)Φ_(θ)(μ,α)],

{circumflex over (Γ)}=Ã⁻¹∂_(z)Φ_(x)(μ,α)[Γ+∂_(z)²Φ_(θ)(μ,α)]∂_(z)Φ_(x)(μ,α)^(T) Ã ^(−T).  (56)

It is noted that (μ,α) is itself a quadrature node corresponding to thecanonical node (z,φ)=(0,0). Hence, {circumflex over (α)} is recovered.The partial derivatives of Φ in the expressions for {circumflex over(β)} and {circumflex over (Γ)} are approximated using a two-pointcentral difference scheme (for the first-order derivatives ∂_(z)Φ_(x)and ∂_(z)Φ_(θ)) and a three-point central difference scheme (for thesecond-order derivatives ∂_(z) ²Φ_(θ)). Shown below are examples using ascalar function ƒ:

${{f^{\prime}(x)} \approx \frac{{f\left( {x + \xi} \right)} - {f\left( {x - \xi} \right)}}{2\xi}},{{f^{''}(x)} \approx {\frac{{f\left( {x + \xi} \right)} - {2{f(x)}} + {f\left( {x - \xi} \right)}}{\xi^{2}}.}}$

The value used for ξ depends on the GVM quadrature rule used to generatethe quadrature nodes. For the third and fifth-order methods, ξ=√{squareroot over (3)}, as given by Equations (39) and (43). No additionalfunction evaluations of Φ are required in these central differencecalculations since they only make use of the transformed quadraturenodes {({tilde over (x)}_(σ) _(i) ,{tilde over (θ)}_(σ) _(i) )}_(i=1)^(N). It is noted that the use of numerical approximations for thepartial derivatives of Φ can be diverted if it is computationallyfeasible to calculate them directly.

When specialized to the two-body problem, some additional assumptionscan be made in this step to reduce computations. If the third-order GVMquadrature method is assumed, then the N=13 canonical sigma points canbe enumerated as

${\begin{bmatrix}z_{\sigma_{1}} & \ldots & z_{\sigma_{13}} \\\varphi_{\sigma_{1}} & \ldots & \varphi_{\sigma_{13}}\end{bmatrix} = \begin{bmatrix}0 & 0 & 0 & \xi & {–\xi} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & \xi & {- \xi} & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & \xi & {- \xi} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \xi & {- \xi} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \xi & {- \xi} \\0 & \eta & {- \eta} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{bmatrix}},$

where ξ and η are specified in (39). Then,

${{\partial_{z}{\Phi_{\theta}\left( {\mu,\alpha} \right)}} \approx {\frac{1}{2\xi}\begin{bmatrix}{{\overset{\sim}{\theta}}_{\sigma_{4}} - {\overset{\sim}{\theta}}_{\sigma_{5}}} \\{{\overset{\sim}{\theta}}_{\sigma_{6}} - {\overset{\sim}{\theta}}_{\sigma_{7}}} \\{{\overset{\sim}{\theta}}_{\sigma_{8}} - {\overset{\sim}{\theta}}_{\sigma_{9}}} \\{{\overset{\sim}{\theta}}_{\sigma_{10}} - {\overset{\sim}{\theta}}_{\sigma_{11}}} \\{{\overset{\sim}{\theta}}_{\sigma_{12}} - {\overset{\sim}{\theta}}_{\sigma_{13}}}\end{bmatrix}}},{{\partial_{z}^{2}{\Phi_{\theta}\left( {\mu,\alpha} \right)}} \approx {{\frac{1}{\xi^{2}}\begin{bmatrix}{{\overset{\sim}{\theta}}_{\sigma_{4}} - {2{\overset{\sim}{\theta}}_{\sigma_{1}}} + {\overset{\sim}{\theta}}_{\sigma_{5}}} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0\end{bmatrix}}.}}$

and Φ_(θ)(μ,α)={tilde over (θ)}_(σ) ₁ . Further, the approximation∂_(x)Φ_(x)≈I is assumed. (These approximations are exact under theassumption of unperturbed two-body dynamics in (5)). Thus,

{circumflex over (α)}=Φ_(θ)(μ,α),{circumflex over (β)}={tilde over(A)}⁻¹ A[β+∂ _(z)Φ_(θ)(μ,α)],{circumflex over (Γ)}={tilde over (A)}⁻¹A[Γ+∂ _(z) ²Φ_(θ)(μ,α)]A ^(T){tilde over (A)}^(−T).

At step 104, the hatted parameters {circumflex over (α)}, {circumflexover (β)}, and {circumflex over (Γ)} computed in Step 103 can beapproximations to the “optimal” parameters {tilde over (α)}, {tilde over(β)}, and {tilde over (Γ)} characterizing the output GVM distribution.In this step, the hatted parameters can be “refined” by solving anon-linear least squares problem which is motivated as follows. Define

p _(a)({tilde over (x)},{tilde over (θ)})=

({tilde over (x)},{tilde over (θ)};{tilde over (μ)},{tilde over(P)},{circumflex over (α)},{circumflex over (β)},{circumflex over(Γ)},{tilde over (κ)}),  (57)

p _(e)({tilde over (x)},{tilde over (θ)})=

(Ψ_(x)({tilde over (x)},{tilde over (θ)};p),Ψ_(θ)({tilde over(x)},{tilde over (θ)};p);μ,P,α,β,Γ,κ.  (58)

It is noted that p_(a) is the approximate GVM distribution of the outputusing the hatted parameters, while p_(e) is the “exact” PDF, as given by(49), but with the assumption that the determinant factor is unity(i.e., Φ is volume preserving). It is acknowledged that this issuboptimal (only if Φ is not volume preserving) though, in the contextof the present invention, it is observed that the largest deviationsfrom unity in the determinant factor are on the order of 10⁻⁴ forscenarios with the strongest non-conservative forces.

In this refinement step, the approximations of the hatted parameters canbe improved by “matching” the approximate and exact PDFs in (57) and(58) at the N transformed quadrature nodes {({tilde over (x)}_(σ) _(i),{tilde over (θ)}_(σ) _(i) )}_(i=1) ^(N). One way to accomplish thisobjective is to study the overdetermined system of algebraic equations

p _(a)({tilde over (x)}_(σ) _(i) ,{tilde over (θ)}_(σ) _(i) )=p_(e)({tilde over (x)}_(σ) _(i) ,{tilde over (θ)}_(σ) _(i) ),i=1, . . .,N,

in {circumflex over (α)}, {circumflex over (β)}, and {circumflex over(Γ)}. Instead, the analogous equations in ‘log space’ are used bydefining

${{l_{a}\left( {\overset{\sim}{x},\overset{\sim}{\theta}} \right)} = {{- 2}{\ln \left\lbrack \frac{p_{a}\left( {\overset{\sim}{x},\overset{\sim}{\theta}} \right)}{p_{a}\left( {{\overset{\sim}{x}}_{*},{\overset{\sim}{\theta}}_{*}} \right)} \right\rbrack}}},{{l_{e}\left( {\overset{\sim}{x},\overset{\sim}{\theta}} \right)} = {{- 2}{\ln \left\lbrack \frac{p_{e}\left( {\overset{\sim}{x},\overset{\sim}{\theta}} \right)}{p_{e}\left( {{\overset{\sim}{x}}_{*},{\overset{\sim}{\theta}}_{*}} \right)} \right\rbrack}}},$

where ({tilde over (x)}_(*),{tilde over (θ)}_(*)) denotes the modalpoint. It follows that

l _(a)({tilde over (x)},{tilde over (θ)})=M({tilde over (x)},{tilde over(θ)};{tilde over (μ)},{tilde over (P)},{circumflex over (α)},{circumflexover (β)},{circumflex over (Γ)},{tilde over (κ)}),

l _(e)({tilde over (x)},{tilde over (θ)})=M(Ψ_(x)({tilde over(x)},{tilde over (θ)};p),Ψ_(θ)({tilde over (x)},{tilde over(θ)};p);μ,P,α,β,Γ,κ),

where M is the Mahalanobis von Mises statistic (26). Now set {tilde over(κ)}=κ (see (53) or (54)) and solve the overdetermined system ofalgebraic equations

l _(a)({tilde over (x)} _(σ) _(i) ,{tilde over (θ)}_(σ) _(i) )=l_(e)({tilde over (x)} _(σ) _(i) ,{tilde over (θ)}_(σ) _(i) ),i=1, . . .,N,

for {circumflex over (α)}, {circumflex over (β)}, and {circumflex over(Γ)}, in the least squares sense. This leads to the optimization problem(55) which can be expressed in the equivalent form

$\begin{matrix}{{{\left( {\overset{\sim}{\alpha},\overset{\sim}{\beta},\overset{\sim}{\Gamma}} \right) = {\begin{matrix}{\arg \mspace{14mu} \min} \\{\hat{\alpha},\hat{\beta},\hat{\Gamma}}\end{matrix}{\sum\limits_{i = 1}^{N}\left\lbrack {r_{i}\left( {\hat{\alpha},\hat{\beta},\hat{\Gamma}} \right)} \right\rbrack^{2}}}},{where}}{{{r_{i}\left( {\hat{\alpha},\hat{\beta},\hat{\Gamma}} \right)} = {{z_{\sigma_{i}}^{T}z_{\sigma_{i}}} - {{\overset{\sim}{z}}_{\sigma_{i}}^{T}{\overset{\sim}{z}}_{\sigma_{i}}} + {4{\kappa \left( {{\sin^{2}\frac{1}{2}\varphi_{\sigma_{i}}} - {\sin^{2}\frac{1}{2}{\overset{\sim}{\varphi}}_{\sigma_{i}}}} \right)}}}},{and}}{{{\overset{\sim}{z}}_{\sigma_{i}} = {{\overset{\sim}{A}}^{- 1}\left( {{\overset{\sim}{x}}_{\sigma_{i}} - \overset{\sim}{\mu}} \right)}},{{\overset{\sim}{\varphi}}_{\sigma_{i}} = {{\overset{\sim}{\theta}}_{\sigma_{i}} - \hat{\alpha} - {{\hat{\beta}}^{T}{\overset{\sim}{z}}_{\sigma_{i}}} - {\frac{1}{2}{\overset{\sim}{z}}_{\sigma_{i}}^{T}\hat{\Gamma}{{\overset{\sim}{z}}_{\sigma_{i}}.}}}}}} & (59)\end{matrix}$

Methods for solving non-linear least squares problems, such asGauss-Newton, full Newton, and quasi-Newton updates, along withglobalization methods such as line search and trust region methodsincluding Levenberg-Marquardt, are efficient and mature and will not bediscussed further here. It is noted that the hatted parameters computedin Step 103 can be used as a starting iteration.

For the two-body problem, it is proposed to optimize over the parameters{circumflex over (α)}, {circumflex over (β)}, and the (1,1)-component of{circumflex over (Γ)}. (The largest changes in the initial parametermatrix T to the transformed matrix {tilde over (Γ)} is in the(1,1)-component.) To facilitate the solution of the non-linear leastsquares problem (59) with these specializations, it is useful to notethe partial derivatives of the residuals r_(i):

${\frac{\partial r_{i}}{\partial\hat{\alpha}} = {2\kappa \; \sin {\overset{\sim}{\varphi}}_{\sigma_{i}}}},{\frac{\partial r_{i}}{\partial\hat{\beta}} = {2{\kappa sin}{\overset{\sim}{\varphi}}_{\sigma_{i}}{\overset{\sim}{z}}_{\sigma_{i}}}},{\frac{\partial r_{i}}{\partial{\hat{\Gamma}}_{11}} = {\kappa \; \sin {{{\overset{\sim}{\varphi}}_{\sigma_{i}}\left\lbrack \left( {\overset{\sim}{z}}_{\sigma_{i}} \right)_{11} \right\rbrack}^{2}.}}}$

It is noted that the optimization problem (59) is numericallywell-conditioned and, for the case of the two-body problem, isinexpensive to solve relative to the propagation of the quadraturenodes.

V.D. Performance Metrics

In this subsection, two different metrics are proposed which allow theoperator to assess online when the non-linear transformation of a GVMdistribution is not well-represented by a GVM distribution. The firstmetric, given by Equation (61), is based on the differential entropy andwas popularized by DeMars et al. in their design of an adaptive Gaussianmixture filter. While the entropy has a nice information-theoreticinterpretation and is straightforward to apply, it is only applicable totransformations which are volume preserving. In the context of thetwo-body problem in orbital mechanics, volume preserving transformationscan be induced generally in the case of pure (conservative) gravityflows (solar radiation pressure is one exception). The second metric,given by Equation (62), is valid for general (non-volume preserving)non-linear transformations and can be based on the Mahalanobis von Misesstatistic. Both metrics provide a tool to validate the GVM propagationalgorithm in the previous subsection. If a breakdown is found to occur,a higher-fidelity method can be triggered such as one that uses amixture of GVM distributions to more accurately characterize thetransformed PDF. Additional details are provided in Subsection V.F.

V.D.i. Differential Entropy

Let x denote a random vector defined on a manifold

with corresponding PDF p_(x)(x). The differential entropy of x isdefined by

$\begin{matrix}{{H(x)} = {- {\int_{\mathcal{M}}{{p_{x}(x)}\ln \; {p_{x}(x)}{{x}.}}}}} & (60)\end{matrix}$

If Φ:

→

is a diffeomorphism that is also volume preserving, i.e.,

${{\det \left\lbrack \frac{\partial\Phi}{\partial x} \right\rbrack} = 1},$

for all xε

, then it follows from the change of variables theorem for PDFs thatH(x)=H({tilde over (x)}), where {tilde over (x)}=Φ(x). In other words,the differential entropy is invariant under a volume preservingdiffeomorphism.

If (x,θ)˜GVM(μ, P, α, β, Γ, κ), then its differential entropy can begiven by the expression in (17). This expression can be used to validatethe approximation ({tilde over (x)},{tilde over (θ)}){dot over(˜)}GVM({tilde over (μ)}, {tilde over (P)}, {tilde over (α)}, {tildeover (β)}, {tilde over (Γ)}, {tilde over (κ)}), where ({tilde over(x)},{tilde over (θ)})=Φ(x,θ) and Φ is a volume preservingdiffeomorphism on

^(n)×

. Indeed, any deviation between the differential entropy of (x,θ) and({tilde over (x)},{tilde over (θ)}) indicates a breakdown of the GVMassumption. To detect such deviations, one can evaluate the percentagechange in differential entropy:

$\begin{matrix}{{\Delta \; {H\left( {\overset{\sim}{x},{\overset{\sim}{\theta};x},\theta} \right)}} = {100{{\frac{{H\left( {\overset{\sim}{x},\overset{\sim}{\theta}} \right)} - {H\left( {x,\theta} \right)}}{H\left( {x,\theta} \right)}}.}}} & (61)\end{matrix}$

If ΔH exceeds a user-defined threshold, then a breakdown in the GVMassumption is declared.

V.D.ii. Mahalanobis von Mises Statistic

Diffeomophisms are generally not volume preserving including thoseencountered in the two-body problem induced by non-conservative forcessuch as atmospheric drag. In such cases, to detect a breakdown in theGVM assumption discussed above and to thereby validate the GVMpropagation algorithm in Subsection V.C, a metric is proposed based onthe Mahalanobis von Mises statistic (26) formed from the quadraturenodes of a GVM quadrature rule.

The metric is motivated from the property that the Mahalanobis von Misesstatistic can be invariant under any transformation of the form (18).Indeed, if (x,θ)˜GVM(μ, P, α, β, Γ, κ), and ({tilde over (x)},{tildeover (θ)}) is defined by the transformation (18), then ({tilde over(x)},{tilde over (θ)})˜GVM({tilde over (μ)}, {tilde over (P)}, {tildeover (α)}, {tilde over (β)}, {tilde over (Γ)}, {tilde over (κ)}), wherethe transformed parameters are given by (19). It follows that

M(x,θ;μ,P,α,β,Γ,κ)=M({tilde over (x)},{tilde over (θ)};{tilde over(μ)},{tilde over (P)},{tilde over (α)},{tilde over (β)},{tilde over(Γ)},{tilde over (κ)}),

where M is the Mahalanobis von Mises statistic (26).

Using the same input from the algorithm of Subsection V.C, define

${\Sigma_{M} = {\sum\limits_{i = 1}^{N}{M\left( {x_{\sigma_{i}},{\theta_{\sigma_{i}};\mu},P,\alpha,\beta,\Gamma,\kappa} \right)}}},$

where {(x_(σ) _(i) ,θ_(σ) _(i) )}_(i=1) ^(N) are the N quadrature nodesgenerated from Step 101 of the algorithm. Equivalently, the expressionfor Σ_(M) can be written in terms of the canonical GVM sigma points{(z_(σ) _(i) ,φ_(σ) _(i) )}_(i=1) ^(N) according to

$\Sigma_{M} = {\sum\limits_{i = 1}^{N}{{M\left( {z_{\sigma_{i}},{\varphi_{i};0},I,0,0,0,\kappa} \right)}.}}$

Next, define

$\begin{matrix}{{{\overset{\sim}{\Sigma}}_{M} = {\sum\limits_{i = 1}^{N}{M\left( {{\overset{\sim}{x}}_{\sigma_{i}},{{\overset{\sim}{\theta}}_{\sigma_{i}};\overset{\sim}{\mu}},\overset{\sim}{P},\overset{\sim}{\alpha},\overset{\sim}{\beta},\overset{\sim}{\Gamma},\overset{\sim}{\kappa}} \right)}}},} & (62)\end{matrix}$

where {({tilde over (x)}_(σ) _(i) ,{tilde over (θ)}_(σ) _(i) )}_(i=1)^(N) are the transformed quadrature nodes computed from Step 101 of thealgorithm and ({tilde over (μ)}, {tilde over (P)}, {tilde over (α)},{tilde over (β)}, {tilde over (Γ)}, {tilde over (κ)}) are the parameterset of the output GVM distribution computed in Steps 102-104.

By the invariance property of the Mahalanobis von Mises statisticdiscussed above, if {tilde over (Σ)}_(M)≠Σ_(M), then the underlyingtransformation that produced the sigma points {({tilde over (x)}_(σ)_(i) ,{tilde over (θ)}_(σ) _(i) )}_(i=1) ^(N) is not of the form (18).Hence, the PDF of the transformed random vector Φ(x,θ) is not the outputGVM distribution of the algorithm. Significant departures of thetransformed PDF from a GVM distribution are manifested in deviations in{tilde over (Σ)}_(M) from Σ_(M). Therefore, if |{tilde over(Σ)}_(M)−Σ_(M)|/Σ_(M) exceeds a user-defined threshold, then a breakdownin the GVM assumption is declared.

V.E. Inclusion of Stochastic Process Noise

The inclusion of uncertain parameters or stochastic process noise in thetransformation of a GVM random vector can be treated within thealgorithm of Subsection V.C using state augmentation. Indeed, let Φ:

^(n)×

^(m)×

→

^(n)×

be a family of diffeomorphisms in wε

^(m) defined such that

({tilde over (x)},{tilde over (θ)})=Φ(x,w,θ;p),

where p denotes any constant non-stochastic parameters. It is assumedthat (y,θ)˜GVM(μ, P, α, β, Γ, κ), where y=(x,w) denotes the augmentedrandom vector. In many filtering applications, it commonly assumed thatthe state vector (x,θ) is independent of the process noise random vectorw and that w is additive so that Φ enjoys the form

Φ(x,w,θ;p)=f(x,θ;p)+g(x,θ;p)w.

In this formulation, such assumptions on the additivity of the processnoise and the independence of (x,θ) with w can be relaxed.

Introducing a slack variable {tilde over (w)}, consider the augmentedtransformation

({tilde over (y)},{tilde over (θ)})=Ω(y,θ;p),  (63)

where {tilde over (y)}=({tilde over (x)},{tilde over (w)}) and

${\Omega \left( {y,\theta} \right)} = {\begin{bmatrix}{\Phi_{x}\left( {x,w,{\theta;p}} \right)} \\w \\{\Phi_{0}\left( {x,w,{\theta;p}} \right)}\end{bmatrix}.}$

It is noted that Ω is a diffeomorphism from

^(n+m)×

to

^(n+m)×

. Therefore, the algorithm of Subsection V.C can be applied inconjunction with the transformation (63) to yield an approximation ofthe joint PDF of ({tilde over (y)},{tilde over (θ)}) by a GVMdistribution. To obtain the joint PDF of ({tilde over (x)},{tilde over(θ)}), the slack variable {tilde over (w)} is marginalized out, notingfrom Section III.G that ({tilde over (x)},{tilde over (θ)}) is also aGVM distribution given that the output ({tilde over (y)},{tilde over(θ)}) is a GVM distribution.

V.F. Propagation of Mixtures

According to one embodiment, the framework of this section can beapplied equally to the propagaAccording to another embodiment, the“adaptive entropy-based Gaussian-mixture information synthesis” (AEGIS)methodology of DeMars et al, a Gaussian mixture filter in which thenumber of mixture components adapts according to the non-linearity, canbe extended to mixtures of GVM distributions using the metrics definedin Subsection V.D.

The random variables (x,θ)ε

^(n)×

are said be jointly distributed as a mixture of GVM distributions iftheir joint PDF has the form

$\begin{matrix}{{{p\left( {x,\theta} \right)} = {\sum\limits_{i = 1}^{N}{w_{i}\left( {x,{\theta;\mu_{i}},P_{i},\alpha_{i},\beta_{i},\Gamma_{i},\kappa_{i}} \right)}}},} & (64)\end{matrix}$

where the weights w_(i), i=1, . . . , N, are positive and sum to unity.Under a diffeomorphism ({tilde over (x)},{tilde over (θ)})=Φ(x,θ; p),the joint PDF of ({tilde over (x)},{tilde over (θ)}) is approximated as

$\begin{matrix}{{{p\left( {\overset{\sim}{x},\overset{\sim}{\theta}} \right)} \approx {\sum\limits_{i = 1}^{N}\; {w_{i}\left( {\overset{\sim}{x},{\overset{\sim}{\theta};{\overset{\sim}{\mu}}_{i}},{\overset{\sim}{P}}_{i},{\overset{\sim}{\alpha}}_{i},{\overset{\sim}{\beta}}_{i},{\overset{\sim}{\Gamma}}_{i},{\overset{\sim}{\kappa}}_{i}} \right)}}},} & (65)\end{matrix}$

where the weights w_(i) can be left unchanged and the parameter sets({tilde over (μ)}_(i), {tilde over (P)}_(i), {tilde over (α)}_(i),{tilde over (β)}_(i), {tilde over (Γ)}_(i), {tilde over (κ)}_(i)) ofeach of the transformed GVM components can be computed using thealgorithm in Subsection V.C. In analogy to the prediction step of theGaussian mixture filter, the procedure is parallelizable since each GVMcomponent can be processed through the algorithm independently of theothers. The approximation in (65) can be provided each input GVMcomponent gets mapped to a GVM component under the transformation Φ. Thefidelity of the approximation can be validated using the metrics inSubsection V.D.

A feature of the Gaussian mixture filter is an algorithm for choosingthe component means, covariances, and weights of the initial inputGaussian mixture so that, when acted on by a non-linear transformation,the resulting output Gaussian mixture is a proper characterization ofthe actual transformed PDF up to a prescribed accuracy. Noting that any(smooth) non-linear transformation is approximately linear in asufficiently small neighborhood, Gaussians with smaller covariancesremain more Gaussian than those with larger covariances under thenon-linear map. Thus, a Gaussian mixture which is refined byapproximating each component Gaussian by a finer Gaussian mixtureexhibits better behavior under a non-linear transformation. Within thetwo-body problem, the refinement is usually performed along onecoordinate direction, namely along the semi-major axis, because it isalong this coordinate where the non-linearities are most severe.

The refinement paradigm for Gaussians and Gaussian mixtures describedabove applies equally well to GVM distributions. Suppose it is desiredto refine the GVM component

${{p\left( {x,\theta} \right)} = {{\left( {x,{\theta;\overset{\_}{\mu}},\overset{\_}{P},\overset{\_}{\alpha},\overset{\_}{\beta},\overset{\_}{\Gamma},\overset{\_}{\kappa}} \right)} = {\left( {{x;\overset{\_}{\mu}},\overset{\_}{P}} \right)\left( {{\theta;{\alpha + {{\overset{\_}{\beta}}^{T}z} + {\frac{1}{2}z^{T}\overset{\_}{\Gamma}\; z}}},\overset{\_}{\kappa}} \right)}}},$

where z=Ā⁻¹(x− μ). If the Gaussian component is refined into a Gaussianmixture according to

${{\left( {{x;\overset{\_}{\mu}},\overset{\_}{P}} \right)} \approx {\sum\limits_{i = 1}^{N}\; {w_{i}\left( {x,\mu_{i},P_{i}} \right)}}},$

then it follows that p(x,θ) is approximated as a mixture of GVMdistributions of the form (64) where

${\alpha_{i} = {\overset{\_}{\alpha} + {{\overset{\_}{\beta}}^{T}v_{i}} + {\frac{1}{2}v_{i}^{T}\overset{\_}{\Gamma}v_{i}}}},{\beta_{i} = {A_{i}^{T}{{\overset{\_}{A}}^{- T}\left( {\overset{\_}{\beta} + {\overset{\_}{\Gamma}v_{i}}} \right)}}},{\Gamma_{i} = {A_{i}^{T}{\overset{\_}{A}}^{- T}\overset{\_}{\Gamma}{\overset{\_}{A}}^{- 1}A_{i}}},{\kappa_{i} = \overset{\_}{\kappa}},{v_{i} = {{\overset{\_}{A}}^{- 1}\left( {\mu_{i} - \overset{\_}{\mu}} \right)}},{{{for}\mspace{14mu} i} = 1},\ldots \mspace{14mu},{N.}$

A feature of the AEGIS methodology for Gaussian mixtures is its abilityto detect online when a single Gaussian or Gaussian mixture, afterundergone a non-linear transformation, does not properly characterizethe actual uncertainty. AEGIS uses the differential entropy (60) todetect departures from “Gaussianity” due to local non-linearities in thetransformation. If a sufficiently large departure is detected, itrefines the Gaussian mixture into a finer Gaussian mixture in order tomitigate the non-linear effects and improve subsequent accuracy. Thus,to adapt the AEGIS concept to mixtures of GVM distributions, thedifferential entropy metric (61), valid for volume preservationtransformations, or the Mahalanobis von Mises metric (62), applicable togeneral non-volume preserving transformations, can be used to detectwhen additional refinement of the input mixture is required.

Stated another way, predicting a location of an object in amulti-dimensional space having a plurality of (n+1) dimensions cancomprise reading 100 a prior probability density function defined by aset of parameters (μ, P, α, β, Γ, κ) and representing an uncertainty ofthe location of the object in a cylindrical manifold (

^(n)×

) of the multi-dimensional space and a diffeomorphism (Φ) of thecylindrical manifold (Φ:

^(n)×

→

^(n)×

). A transformed probability density function defined by a set ofparameters ({tilde over (μ)}, {tilde over (P)}, {tilde over (α)}, {tildeover (β)}, {tilde over (Γ)}, {tilde over (κ)}) and representing anuncertainty of the location of the object in a cylindrical manifold canbe generated 101-104. The transformed probability density function canbe generated 101-104 from the input probability density function underthe diffeomorphism. The input probability density function and thetransformed probability density function are both represented by Gaussvon Mises distributions defined on the cylindrical manifold. The set ofparameters representing the transformed probability density function canthen be provided 105 as a representation of the predicted location ofthe object in the multi-dimensional space. In a SSA application, thediffeomorphism can be a solution flow induced from a system of ordinarydifferential equations describing two-body dynamics of orbital mechanicsfor the object. Also in such an application, the input Gauss von Misesdistribution can be generated from a plurality of radar,electro-optical, or infrared sensor observations.

Generating the transformed probability density function can comprisecomputing 102 the parameters {tilde over (μ)} and {tilde over (P)} froma sequence of function evaluations Φ(x_(σ) _(i) ,θ_(σ) _(i) ), for i=1,. . . , N, where (x_(σ) _(i) ,θ_(σ) _(i) ), for i=1, . . . , N, are achosen sequence of quadrature nodes on

^(n)×

and wherein the quadrature nodes (x_(σ) _(i) , θ_(σ) _(i) ), for i=1, .. . , N, are generated 101 from a Gauss von Mises quadrature rule of achosen order of accuracy in conjunction with the input parameter set (μ,P, α, β, Γ, κ). Additionally, {tilde over (κ)} can be set to {tilde over(κ)}=κ. The parameters {tilde over (α)}, {tilde over (β)}, and {tildeover (Γ)} by {circumflex over (α)}, {circumflex over (β)}, and{circumflex over (Γ)}, respectively, can be approximated 103 usingexpressions depending on the partial derivatives of the diffeomorphismΦ. The parameters {tilde over (α)}, {tilde over (β)}, and {tilde over(Γ)} can be selected 104 according to:

${\left( {\overset{\sim}{\alpha},\overset{\sim}{\beta},\overset{\sim}{\Gamma}} \right) = {\begin{matrix}{\arg {\; \;}\min} \\{\hat{\alpha},\hat{\beta},\hat{\Gamma}}\end{matrix}{\sum\limits_{i = 1}^{N}\; \left\lbrack {{M\left( {x_{\sigma_{i}},{\theta_{\sigma_{i}};\mu},P,\alpha,\beta,\Gamma,\kappa} \right)} - {M\left( {{{\Phi \left( {x_{\sigma_{i}},\theta_{\sigma_{i}}} \right)};\overset{\sim}{\mu}},\overset{\sim}{P},\hat{\alpha},\hat{\beta},\hat{\Gamma},\overset{\sim}{\kappa}} \right)}} \right\rbrack^{2}}}},$

where M is a Mahalanobis von Mises statistic. In some cases, theaccuracy of the transformed Gauss von Mises distribution can bevalidated using a differential entropy metric. Additionally oralternatively, the accuracy of the transformed Gauss von Misesdistribution can be validated using a metric based on a Mahalanobis vonMises statistic.

VI. Data Fusion

The problem of data fusion in multi-target tracking entails thecombining of reports emanating from a common object in order to improvethe state or understanding of that object. The “combining” or “fusion”process is performed in a probabilistic manner based on the laws ofconditional probability and Bayes' rule. Within a sequential non-linearfilter, this operation is called the correction or measurement updatestep of the filter.

This section shows how the correction step of the Bayesian non-linearfilter can be specialized to the case when the prior state is a randomvector represented by a Gauss von Mises (GVM) distribution and theupdate report, hypothesized to emanate from the prior, is either (i)another GVM random vector of the same dimension as the prior, or (ii) anobservation related to the prior by a stochastic measurement model. Ineither case, the output of the Bayesian correction step can be a GVMdistribution characterizing the fusion of the prior and the update. Inaddition, the process can furnish a statistically rigorous predictionerror, a term appearing in the likelihood ratios for scoring theassociation of one report to another. The algorithm descriptions forCases (i) and (ii) are provided in Subsections VI.A and VI.B,respectively, with respective block diagrams provided in FIGS. 2 and 3.

VI.A. Full State Update

FIG. 2 is a flowchart illustrating an exemplary process for fusing aprior state represented by a Gauss von Mises distribution with an updatereport, wherein said update is another Gauss von Mises distribution ofthe same dimension as the prior according to one embodiment of thepresent invention. In the case of a full state update, the Bayesiannon-linear filter correction step has two inputs 200:

-   -   i. a random vector (x₁,θ₁)ε        ^(n)×        , called the prior, distributed as a GVM distribution with        parameter set (μ₁, P₁, α₁, β₁, Γ₁, κ₁), and    -   ii. a second random vector (x₂,θ₂)ε        ^(n)×        , called the update, distributed as a GVM distribution with        parameter set (μ₂, P₂, α₂, β₂, Γ₂, κ₂).        There are two outputs 203:    -   i. a random vector (x,θ)ε        ^(n)×        , called the correction, distributed as a GVM distribution with        parameter set (μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ)) which        characterizes the fusion of the prior with the update, and    -   ii. the prediction error used to score the association of the        prior with the update.        The algorithm can have two main steps summarized below:    -   Compute 201 the parameter set (μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ),        Γ_(ƒ), κ_(ƒ)) of the corrected GVM distribution from        Equations (71) and (74); and    -   Evaluate 202 the prediction error integral (75) using a GVM        quadrature rule in conjunction with the output from Step 201 and        Equation (76).        These two steps 201 and 202 are described in more detail below        including discussions on theoretical considerations.

VI.A.i. Theoretical Considerations

In a general setting, suppose X and Y are random vectors defined on amanifold M with respective probability density functions (PDFs) p_(X|I)_(x) (x|I_(x)) and p_(Y|I) _(y) (y|I_(y)), where I_(x) and I_(y) denoteany prior information associated with X and Y, respectively. Denote thejoint PDF of X and Y (conditioned on their respective prior information)by p_(X,Y|I) _(x) _(,I) _(y) (x, y|I_(x), I_(y)) and the conditional PDFof X given Y=y, I_(x), and I_(y) by p_(X|Y,I) _(x) _(,I) _(y) (x|y,I_(x), I_(y)).

The PDF obtained by fusing the information of X with that of Y isdefined to be the conditional PDF of X given Y=x and given the otherprior information:

p _(X→Y|I) _(x) _(,I) _(y) (x|I _(x) ,I _(y))≡p _(X|Y,I) _(x) _(,I) _(y)(x|x,I _(x) ,I _(y)).

By the definition of conditional probability,

$\begin{matrix}{{{p_{{{X->Y}I_{x}},I_{y}}\left( {{xI_{x}},I_{y}} \right)} = {\frac{1}{c}{p_{X,{YI_{x}},I_{y}}\left( {x,{xI_{x}},I_{y}} \right)}}},} & (66)\end{matrix}$

where the normalization constant c is given by

c=p _(Y|I) _(x) _(,I) _(y) (x|I _(x) ,I _(y))=

p_(X,Y|I) _(x) _(,I) _(y) (x,x|I _(x) ,I _(y))dx  (67)

Additionally, if (i) X is independent of Y given the prior informationI_(x) and I_(y), (ii) X is independent of I_(y) given I_(x), and (iii) Yis independent of I_(x) given I_(y), then (66) and (67) simplify to

${{p_{{{X->Y}I_{x}},I_{y}}\left( {{xI_{x}},I_{y}} \right)} = {\frac{1}{c}{p_{XI_{x}}\left( {xI_{x}} \right)}{p_{YI_{y}}\left( {xI_{y}} \right)}}},{c = {{p_{XI_{x}}\left( {xI_{x}} \right)}{p_{YI_{y}}\left( {xI_{y}} \right)}\ {{x}.}}}$

In what follows, the above theoretical considerations and independenceassumptions are specialized to the case when the random vectors X and Yare the input described at the beginning of the subsection (and I_(x)and I_(y) denote any prior information used to construct theirrespective parameter sets).

VI.A.ii. Algorithm Description

At step 201, the fused (corrected) PDF in (x,θ) obtained by fusing theprior random vector (x₁,θ₁) with the information from the update randomvector (x²,θ₂) is

$\begin{matrix}{{{p_{f}\left( {x,\theta} \right)} = {\frac{1}{c}\left( {x,{\theta;\mu_{1}},P_{1},\alpha_{1},\beta_{1},\Gamma_{1},\kappa_{1}} \right)\left( {x,{\theta;\mu_{2}},P_{2},\alpha_{2},\beta_{2},\Gamma_{2},\kappa_{2}} \right)}},} & (68) \\{\mspace{79mu} {where}} & \; \\{c = {\int_{{\mathbb{R}}^{n}}{\int_{- \pi}^{\pi}{\ \left( {x,{\theta;\mu_{1}},P_{1},\alpha_{1},\beta_{1},\Gamma_{1},\kappa_{1}} \right)\left( {x,{\theta;\mu_{2}},P_{2},\alpha_{2},\beta_{2},\Gamma_{2},\kappa_{2}} \right){\theta}\ {{x}.}}}}} & (69)\end{matrix}$

The normalization constant c is what is referred to as the predictionerror. The objective is to approximate (68) by a GVM distribution withparameters (μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ)). To proceed, write

$\begin{matrix}{{{{p_{f}\left( {x,\theta} \right)} = {\propto {\left( {{x;\mu_{1}},P_{1}} \right)\left( {{x;\mu_{2}},P_{2}} \right){\exp \left\lbrack {{\kappa_{1}{\cos \left( {\theta - {\Theta_{1}(x)}} \right)}} + {\kappa_{2}{\cos \left( {\theta - {\Theta_{2}(x)}} \right)}}} \right\rbrack}}}},\mspace{79mu} {where}}\mspace{79mu} {{{\Theta_{1}(x)} = {\alpha_{1} + {\beta_{1}^{T}z_{1}} + {\frac{1}{2}z_{1}^{T}\Gamma_{1}z_{1}}}},{z_{1} = {A_{1}^{- 1}\left( {x - \mu_{1}} \right)}},\mspace{79mu} {{\Theta_{2}(x)} = {\alpha_{2} + {\beta_{2}^{T}z_{2}} + {\frac{1}{2}z_{2}^{T}\Gamma_{2}z_{2}}}},{z_{2} = {A_{2}^{- 1}\left( {x - \mu_{2}} \right)}},}} & (70)\end{matrix}$

and ‘∝’ denotes equality up to an overall multiplicative constant.Combining the two Gaussian PDFs in (70) into a single Gaussian PDF, itfollows that

p _(ƒ)(x,θ)∝

(x;μ _(ƒ) ,P _(ƒ))exp[κ₁ cos(θ−Θ(x))+κ₂ cos(θ−Θ₂(x))],

where μ_(ƒ) and P_(ƒ) satisfy

P _(ƒ) ⁻¹ =P ₁ ⁻¹ +P ₂ ⁻¹ ,P _(ƒ) ⁻¹μ_(ƒ) =P ₁ ⁻¹μ₁ +P ₂ ⁻¹μ₂.  (71)

Using elementary trigonometric identities, it follows that

     κ₁cos (θ − Θ₁(x)) + κ₂cos (θ − Θ₂(x)) = K(x)cos (θ − Q(x)),      where$\mspace{79mu} {{{K(x)} = \left\lbrack {\kappa_{1}^{2} + {2\kappa_{1}\kappa_{2}{\cos \left( {{\Theta_{1}(x)} - {\Theta_{2}(x)}} \right)}} + \kappa_{2}^{2}} \right\rbrack^{\frac{1}{2}}},{{Q(x)} = {{\arctan \left( {{{\kappa_{1}\sin \; {\Theta_{1}(x)}} + {\kappa_{2}\sin \; {\Theta_{2}(x)}}},{{\kappa_{1}\cos \; {\Theta_{1}(x)}} + {\kappa_{2}\cos \; {\Theta_{2}(x)}}}} \right)}.}}}$

(The two-argument inverse tangent function is implied.) The final stepin the computation is to find the parameters (α_(ƒ), β_(ƒ), Γ_(ƒ),κ_(ƒ)) such that

K(x)cos(θ−Q(x))≈κ_(ƒ)cos(θ−Θ_(ƒ)(x)),  (72)

where Θ_(ƒ)(x)=α_(ƒ)+β_(ƒ) ^(T)z_(ƒ)+1/2z_(ƒ) ^(T)Γ_(ƒ)z_(ƒ),z_(ƒ)=A_(ƒ) ⁻¹(x−μ_(ƒ)), and P_(ƒ)=A_(ƒ)A_(ƒ) ^(T). To streamlinesubsequent notation, define the auxiliary quantities

$\begin{matrix}{{v_{i} = {A_{i}^{- 1}\left( {\mu_{f} - \mu_{i}} \right)}},{\theta_{i} = {{\Theta_{i}\left( \mu_{f} \right)} = {\alpha_{i} + {\beta_{i}^{T}v_{i}} + {\frac{1}{2}v_{i}^{T}\Gamma_{i}v_{i}}}}},{b_{i} = {A_{i}^{- T}\left( {\beta_{i} + {\Gamma_{i}v_{i}}} \right)}},{C_{i} = {A_{i}^{- T}\Gamma_{i}A_{i}^{- 1}}},} & (73)\end{matrix}$

for i=1,2. The approximation (72) can be achieved using Taylorexpansions yielding

$\begin{matrix}{\mspace{79mu} {{{\kappa_{f} = {{K\left( \mu_{f} \right)} = \left\lbrack {\kappa_{1}^{2} + {2\kappa_{1}\kappa_{2}{\cos \left( {\theta_{1} - \theta_{2}} \right)}} + \kappa_{2}^{2}} \right\rbrack^{1/2}}},\begin{matrix}{\mspace{79mu} {\alpha_{f} = {Q\left( \mu_{f} \right)}}} \\{{= {\arctan \left( {{{\kappa_{1}\sin \; \theta_{1}} + {\kappa_{2}\sin \; \theta_{2}}},{{\kappa_{1}\cos \; \theta_{1}} + {\kappa_{2}\cos \; \theta_{2}}}} \right)}},}\end{matrix}}\mspace{79mu} {{\beta_{f} = {A_{f}^{T}b_{f}}},\mspace{79mu} {\Gamma_{f} = {A_{f}^{T}C_{f}A_{f}}},\mspace{79mu} {where}}{{b_{f} = {\left\lbrack {{\kappa_{1}{\cos \left( {\alpha_{f} - \theta_{1}} \right)}b_{1}} + {\kappa_{2}{\cos \left( {\alpha_{f} - \theta_{2}} \right)}b_{2}}} \right\rbrack/\left\lbrack {{\kappa_{1}{\cos \left( {\alpha_{f} - \theta_{1}} \right)}} + {\kappa_{2}{\cos \left( {\alpha_{f} - \theta_{2}} \right)}}} \right\rbrack}},{C_{f} = {\left\{ {{\kappa_{1}{{\sin \left( {\alpha_{f} - \theta_{1}} \right)}\left\lbrack {{b_{1}b_{1}^{T}} - \left( {{b_{1}b_{f}^{T}} + {b_{f}b_{1}^{T}}} \right) + {b_{f}b_{f}^{T}}} \right\rbrack}} + {\kappa_{2}{{\sin \left( {\alpha_{f} - \theta_{2}} \right)}\left\lbrack {{b_{2}b_{2}^{T}} - \left( {{b_{2}b_{f}^{T}} + {b_{f}b_{2}^{T}}} \right) + {b_{f}b_{f}^{T}}} \right\rbrack}} + {\kappa_{1}{\cos \left( {\alpha_{f} - \theta_{1}} \right)}C_{1}} + {\kappa_{2}{\cos \left( {\alpha_{f} - \theta_{2}} \right)}C_{2}}} \right\}/{\left\lbrack {{\kappa_{1}{\cos \left( {\alpha_{f} - \theta_{1}} \right)}} + {\kappa_{2}{\cos \left( {\alpha_{f} - \theta_{2}} \right)}}} \right\rbrack.}}}}}} & (74)\end{matrix}$

In summary, Equations (71) and (74) specify the parameter set (μ_(ƒ),P_(ƒ), α_(ƒ), β^(ƒ), Γ_(ƒ), κ_(ƒ)) of the output GVM distribution.

At step 202, using the GVM distribution with parameter set (μ_(ƒ),P_(ƒ), α_(ƒ), β^(ƒ), Γ_(ƒ), κ_(ƒ)) computed from above, the predictionerror c given by (69) can be expressed in the form

$\begin{matrix}{{{c = {k{\int_{{\mathbb{R}}^{n}}{\int_{- \pi}^{\pi}{\ \left( {x,{\theta;\mu_{f}},P_{f},\alpha_{f},\beta_{f},\Gamma_{f},\kappa_{f}} \right){f\left( {x,\theta} \right)}{\theta}\ {x}}}}}},\mspace{79mu} {where}}\mspace{79mu} {k = {\frac{2{\pi }^{{- \kappa}\; f}{I_{0}\left( \kappa_{f} \right)}}{2{\pi }^{- {\kappa \;}_{1}}{{I_{0}\left( \kappa_{1} \right)} \cdot 2}{\pi }^{- {\kappa \;}_{2}}{I_{0}\left( \kappa_{2} \right)}} \cdot \left\lbrack \frac{\det \left( {2\pi \; P_{f}} \right)}{{\det \left( {2\pi \; P_{1}} \right)}{\det \left( {2\pi \; P_{2}} \right)}} \right\rbrack^{1/2}}}\mspace{79mu} {and}{{f\left( {x,\theta} \right)} = {{\exp\left\lbrack {{{- \frac{1}{2}}z_{1}^{T}z_{1}} - {\frac{1}{2}z_{2}^{T}z_{2}} + {\frac{1}{2}z_{f}^{T}z_{f}} - {2\kappa_{1}\sin^{2}\frac{1}{2}\left( {\theta - {\Theta_{1}(x)}} \right)} - {2\kappa_{2}\sin^{2}\frac{1}{2}\left( {\theta - {\Theta_{2}(x)}} \right)} + {2\kappa_{f}\sin^{2}\frac{1}{2}\left( {\theta - {\Theta_{f}(x)}} \right)}} \right\rbrack}.}}} & (75)\end{matrix}$

The definitions of z₁, z₂, z_(ƒ), Θ₁(x), Θ₂(x), and Θ_(ƒ)(x) areprovided following Equations (70) and (72). The framework of Gauss vonMises quadrature, as developed in Section IV, is applicable to theevaluation of the integral in (75). If {(z_(σ) _(i) , φ_(σ) _(i))}_(i=1) ^(N) and {w_(σ) _(i) }_(i=1) ^(N) are the quadrature nodes(sigma points) and weights for the canonical GVM distribution, then itfollows that

$\begin{matrix}{\mspace{79mu} {{{c \approx {k{\sum\limits_{i = 1}^{N}\; {w_{\sigma_{i}}f_{\sigma_{i}}}}}},\mspace{79mu} {where}}{{f_{\sigma_{i}} = {\exp\left\lbrack {{{- \frac{1}{2}}v_{1,\sigma_{i}}^{T}v_{1,\sigma_{i}}} - {\frac{1}{2}v_{2,\sigma_{i}}^{T}v_{2,\sigma_{i}}} + {\frac{1}{2}z_{\sigma_{i}}^{T}z_{\sigma_{i}}} - {2\kappa_{1}\sin^{2}\frac{1}{2}\eta_{1,\sigma_{i}}} - {2\kappa_{2}\sin^{2}\frac{1}{2}\eta_{2,\sigma_{i}}} + {2\kappa_{f}\sin^{2}\frac{1}{2}\varphi_{\sigma_{i}}}} \right\rbrack}},\mspace{79mu} {and}}\mspace{79mu} {{v_{j,\sigma_{i}} = {A_{j}^{- 1}\left( {\mu_{f} - \mu_{j} + {A_{f}z_{\sigma_{i}}}} \right)}},{\eta_{j,\sigma_{i}} = {\left( {\varphi_{\sigma_{i}} + \alpha_{f} + {\beta_{f}^{T}z_{\sigma_{i}}} + {\frac{1}{2}z_{\sigma_{i}}^{T}\Gamma_{f}z_{\sigma_{i}}}} \right) - \left( {\alpha_{j} + {\beta_{j}^{T}v_{j,\sigma_{i}}} + {\frac{1}{2}v_{j,\sigma_{i}}^{T}\Gamma_{j}v_{j,\sigma_{i}}}} \right)}},{\quad\mspace{79mu} {{{{for}\mspace{14mu} j} = 1},2.}}}}} & (76)\end{matrix}$

It is noted that the evaluation of the prediction error integral (75)using different order GVM quadrature rules allows the operator to assesswhen the fused PDF (68) is not well-represented by a GVM distributionwith the parameter set (μ_(ƒ), P_(ƒ), α_(ƒ), β^(ƒ), Γ_(ƒ), κ_(ƒ)).Indeed, the approximation (76) to the prediction error is exact if thefused PDF (68) is identically a GVM distribution. In such a case, thefunction ƒ(x,θ) in (75) is a constant. Hence, the prediction errorintegral can be evaluated exactly using any GVM quadrature rule. Forexample, using a trivial first-order GVM quadrature rule leads to

${c_{1} = {k\mspace{11mu} {\exp\left\lbrack {{{- \frac{1}{2}}v_{1}^{T}v_{1}} - {\frac{1}{2}v_{2}^{T}v_{2}} - {2\kappa_{1}\sin^{2}\frac{1}{2}\left( {\alpha_{f} - \theta_{1}} \right)} - {2\kappa_{2}\sin^{2}\frac{1}{2}\left( {\alpha_{f} - \theta_{2}} \right)}} \right\rbrack}}},$

where ν₁, ν₂, θ₁, and θ₂ are the auxiliary quantities defined in (73).(The subscript ‘1’ underneath c implies the use of a first-order GVMquadrature rule.) Similarly, one can evaluate (75) using a third-orderquadrature method in conjunction with Equations (76). The (absolute orrelative) difference between the high and low-order approximations tothe prediction error (75) can be used to estimate the error of theintegration and validate the underlying assumptions of the algorithm. Abreakdown can be declared if this difference exceeds a user-definedthreshold.

VI.B. Measurement-Based Update

FIG. 3 is a flowchart illustrating an exemplary process for fusing aprior state represented by a Gauss von Mises distribution with an updatereport, wherein said update is an observation related to the prior by astochastic measurement model according to one embodiment of the presentinvention. In the case of a measurement-based update, the Bayesiannon-linear filter correction step can have four inputs 300:

-   -   i. a random vector (x,θ) ε        R^(n)×        , called the prior, distributed as a GVM distribution with        parameter set (μ, P, α, β, Γ, κ),    -   ii. a smooth map h:        ^(n)×        →        ^(p), called the measurement model,    -   iii. a random vector vε        ^(p), called the measurement error, distributed as a zero-mean        multivariate Gaussian random vector with covariance R, and    -   iv. a vector yε        ^(p), called the measurement update, hypothesized to be a        realization of the random vector h(x,θ)+v.        There are two outputs 305:    -   i. a random vector (x,θ) ε        ^(n)×        , called the correction, distributed as a GVM distribution with        parameter set (μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ)) which        characterizes the fusion of the prior with the measurement        update, and    -   ii. the prediction error used to score the association of the        prior with the measurement update.        The algorithm can have four main steps summarized below:    -   Compute 301 the components μ_(ƒ) and α_(ƒ) of the output GVM        parameter set by solving the least squares problem (81);    -   Recover 302 the components P_(ƒ), β_(ƒ), and κ_(ƒ) of the output        GVM parameter set from Equation (82);    -   Compute 303 the component Γ_(ƒ) of the output GVM parameter set        by solving the least squares problem (83). If it is believed        that the output is well-approximated by a Gaussian, this step        can be skipped and one can set Γ_(ƒ)=0; and    -   Evaluate 304 the prediction error integral (84) using a GVM        quadrature rule in conjunction with the output from Steps        301-303 and Equation (85).        These four steps are described in more detail below including        discussions on theoretical considerations.

VI.B.i. Theoretical Considerations

In a general setting, suppose X and V are random vectors defined onmanifolds

and

, respectively, with respective PDFs p_(X|I) _(x) (x|I_(x)) andp_(V)(v), where I_(x) denotes any prior information associated with X.Given a smooth map h:

→

, define the random vector Y according to

Y=h(X)+V.

The PDF obtained by fusing the information of X with that of a givenupdate yε

is defined to be the conditional PDF of X given that Y=y and given theother prior information:

p _(X→Y|I) _(x) (x|y,I _(x))≡p _(X|Y,I) _(x) (x|y,I _(x)).

By Bayes' rule and the preceding definitions, it follows that

$\begin{matrix}{{{p_{X->Y}\left( {\left. x \middle| y \right.,I_{x}} \right)} = {\frac{{p_{{Y|X},I_{x}}\left( {\left. y \middle| x \right.,I_{x}} \right)}{p_{X|I_{x}}\left( x \middle| I_{x} \right)}}{p_{Y|I_{x}}\left( y \middle| I_{x} \right)} = {\frac{1}{c}{p_{v}\left( {y - {h(x)}} \right)}{p_{X|I_{x}}\left( x \middle| I_{x} \right)}}}},} & (77)\end{matrix}$

where the normalization c is given by

c=p _(Y|I) _(x) (y|I _(x))=

p_(V)(y−h(x))p _(X|I) _(x) (x|I _(x))dx.  (78)

The above theoretical considerations are now specialized to the casewhen the random vectors X and V and the map h are the input described atthe beginning of the subsection.

At step 301, the fused (corrected) PDF in (x,θ) obtained by fusing theprior GVM distribution with the information from the measurement updatey is

$\begin{matrix}{{{p_{f}\left( {x,\theta} \right)} = {\frac{1}{c}\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right){N\left( {{{y - {h\left( {x,\theta} \right)}};0},R} \right)}}},\mspace{20mu} {where}} & (79) \\{c = {\int_{{\mathbb{R}}^{n}}{\int_{- \pi}^{\pi}{\left( {x,{\theta;\mu},P,\alpha,\beta,\Gamma,\kappa} \right){N\left( {{{y - {h\left( {x,\theta} \right)}};0},R} \right)}{\theta}{{x}.}}}}} & (80)\end{matrix}$

The normalization constant c is what is referred to as the predictionerror. The objective is to approximate (79) by a GVM distribution withparameters (μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ)). The mode of thisoutput GVM distribution is (μ_(ƒ), α_(ƒ)). These parameters can bechosen to coincide with the mode of (79). The computation of the modecan be formulated in terms of a non-linear least squares problem asfollows. Define

${{{q_{f}\left( {x,\theta} \right)} \equiv {- {{\ln p}_{f}\left( {x,\theta} \right)}}} = {{\frac{1}{2}{r\left( {x,\theta} \right)}^{T}{r\left( {x,\theta} \right)}} + k}},$

where k is a constant (independent of x and θ), and

${{r\left( {x,\theta} \right)} = \begin{bmatrix}z \\{2\sqrt{\kappa}\sin \frac{1}{2}\varphi} \\\xi\end{bmatrix}},{where}$${z = {A^{- 1}\left( {x - \mu} \right)}},{\varphi = {\theta - \alpha - {\beta^{T}z} - {\frac{1}{2}z^{T}\Gamma \; z}}},{\xi = {S^{- 1}\left( {{h\left( {x,\theta} \right)} - y} \right)}},{\xi = {S^{- T}\xi}},{R = {{SS}^{T}.}}$

Thus, the parameters μ_(ƒ) and α_(ƒ) can be computed by solving thenon-linear least squares problem

$\begin{matrix}{\left( {\mu_{f},\alpha_{f}} \right) = {{\underset{x,\theta}{argmax}{p_{f}\left( {x,\theta} \right)}} = {\underset{x,\theta}{argmin}\frac{1}{2}{r\left( {x,\theta} \right)}^{T}{{r\left( {x,\theta} \right)}.}}}} & (81)\end{matrix}$

At step 302, the parameters P_(ƒ), P_(ƒ), and K_(ƒ) of the output GVMdistribution can be recovered from the osculating Gaussian of (79) inanalogy to the derivation of Equation (24) in Section III.H.Specifically,

$\begin{matrix}{{\begin{bmatrix}{{A_{f}^{- T}\left( {I + {\kappa_{f}\beta_{f}\beta_{f}^{T}}} \right)}A_{f}^{- 1}} & {{- \kappa_{f}}A_{f}^{- T}\beta_{f}} \\{{- \kappa_{f}}\beta_{f}^{T}A_{f}^{- 1}} & \kappa_{f}\end{bmatrix} = {\partial_{x}^{2}{q_{f}\left( {\mu_{f},\alpha_{f}} \right)}}},} & (82)\end{matrix}$

where P_(ƒ)=A_(ƒ)A_(ƒ) ^(T) and

=(x,θ). Analytic expressions for the second derivatives of q_(ƒ) aregiven by

${{\partial_{x}^{2}{q_{f}\left( {x,\theta} \right)}} = {{{\partial_{x}{r\left( {x,\theta} \right)}^{T}}{\partial_{x}{r\left( {x,\theta} \right)}}} + {\sum\limits_{i = 1}^{n + 1 + p}{{r_{i}\left( {x,\theta} \right)}{\partial_{x}^{2}{r_{i}\left( {x,\theta} \right)}}}}}},$

where r_(i)(x,θ) denotes the i-th component of the residual r(x,θ). TheJacobian of r with respect to

is

${\partial_{x}{r\left( {x,\theta} \right)}} = {\begin{bmatrix}A^{- 1} & 0 \\{{- \sqrt{\kappa}}\cos \frac{1}{2}{\varphi \left( {\beta + {\Gamma \; z}} \right)}^{T}A^{- 1}} & {\sqrt{\kappa}\cos \frac{1}{2}\varphi} \\{S^{- 1}{\partial_{x}{h\left( {x,\theta} \right)}}} & {S^{- 1}{\partial_{\theta}{h\left( {x,\theta} \right)}}}\end{bmatrix}.}$

For the second derivatives, it is first first noted that ∂

²r_(i)(x,θ)=∂

²z_(i)=0, for i=1, . . . , n. For i=π+1,

${\partial_{x}^{2}\left( {2\sqrt{\kappa}\sin \frac{1}{2}\varphi} \right)} = {{\sqrt{\kappa}\begin{bmatrix}{{- {A^{- T}\begin{bmatrix}{{\cos \frac{1}{2}{\varphi\Gamma}} +} \\{\frac{1}{2}\sin \frac{1}{2}{\varphi \left( {\beta + {\Gamma \; z}} \right)}\left( {\beta + {\Gamma \; z}} \right)^{T}}\end{bmatrix}}}A^{- 1}} & {\frac{1}{2}\sin \frac{1}{2}\varphi \; {A^{- 1}\left( {\beta + {\Gamma \; z}} \right)}} \\{\frac{1}{2}\sin \frac{1}{2}{\varphi \left( {\beta + {\Gamma \; z}} \right)}^{T}A^{- 1}} & {{- \frac{1}{2}}\sin \frac{1}{2}\varphi}\end{bmatrix}}.}$

Finally, for r_(i)=ξ_(i),

${\sum\limits_{i = {n + 2}}^{n + 2 + p}{{r_{i}\left( {x,\theta} \right)}{\partial_{x}^{2}{r_{i}\left( {x,\theta} \right)}}}} = {{\sum\limits_{j = 1}^{p}{\xi_{i}{\partial_{x}^{2}\xi_{i}}}} = {\sum\limits_{j = 1}^{p}{{\hat{\xi}}_{j}{{\partial_{x}^{2}{h_{j}\left( {x,\theta} \right)}}.}}}}$

These partial derivatives are also useful if solving the optimizationproblem (81) in Step 301 using Gauss-Newton or full Newton iterations.

At step 303, in the space surveillance tracking problem, the prior GVMdistribution used as an input to this algorithm can be significantlynon-Gaussian characterized by “banana” or “boomerang” shaped level setsand with a parameter matrix Γ possessing non-negligible components. This“non-Gaussianity” in the input is often induced from the propagation ofanother GVM distribution under the non-linear dynamics of the two-bodyproblem. Further, in this application, the measurement update can oftenbe a radar, electro-optical, or infrared sensor observation where thecorresponding measurement model is the map from (equinoctial orbitalelement) state space to the sensor space coordinates. Upon performingthe Bayesian correction step with this input, it is often found that thecorrected (fused) distribution collapses to something which is nearlyGaussian. In other words, ∥Γ_(ƒ)∥≈0. In general, if it is believed thatthe output is well-approximated by a Gaussian, this step can be skippedand one can set Γ_(ƒ)=0. Regardless of whether this step is executed,some insight on how to validate the realism of the corrected GVMdistribution compared to the actual corrected PDF are provided in thediscussions at the end of Step 304.

Let {circumflex over (Γ)}_(ƒ) be an approximation to Γ_(ƒ) and initiallyset {circumflex over (Γ)}_(ƒ)=0. In this step, the approximation{circumflex over (Γ)}_(ƒ) can be refined by solving a non-linear leastsquares problem, in analogy to Step 104 of the uncertainty propagationalgorithm in Section V.C. Using the exact expression for the corrected(fused) PDF p_(ƒ)(x,θ) given by (79) and the approximation given by

p _(a)(x,θ)=

(x,θ;μ _(ƒ) ,P _(ƒ),α_(ƒ),β_(ƒ),{circumflex over (Γ)}_(ƒ),κ_(ƒ)),

the approximation of {circumflex over (Γ)}_(ƒ) is improved by “matching”p_(a) and P_(ƒ) at the N quadrature nodes {(x_(σ) _(i) , θ_(σ) _(i)}_(i=1) ^(N) of the GVM distribution with parameter set (μ_(ƒ), P_(ƒ),α_(ƒ), β_(ƒ), {circumflex over (Γ)}_(ƒ), κ_(ƒ)). These sigma points aredefined from either the third, fifth, or higher-order GVM quadraturerules derived in Sections IV.A, IV.B, and IV.C, respectively. It isnoted that

${x_{\sigma_{i}} = {\mu_{f} + {A_{f}z_{\sigma_{i}}}}},{\theta_{\sigma_{i}} = {\varphi_{\sigma_{i}} + \alpha_{f} + {\beta_{f}^{T}z_{\sigma_{i}}} + {\frac{1}{2}z_{\sigma_{i}}^{T}{\hat{\Gamma}}_{f}z_{\sigma_{i}}}}},$

for i=1, . . . , N, where {(z_(σ) _(i) ,φ_(σ) _(i) )}_(i=1) ^(N) are thequadrature nodes corresponding to the canonical GVM distribution. Oneway to perform this matching procedure is to study the overdeterminedsystem of algebraic equations

p _(ƒ)(x _(σ) _(i) ,θ_(σ) _(i) )=p _(a)(x _(σ) _(i) ,θ_(σ) _(i) ),i=1, .. . ,N,

in {circumflex over (Γ)}_(ƒ). Instead, the analogous equations in ‘logspace’ are used by defining

${{_{f}\left( {x,\theta} \right)} = {{- 2}{\ln \left\lbrack \frac{p_{f}\left( {x,\theta} \right)}{p_{f}\left( {\mu_{f},\alpha_{f}} \right)} \right\rbrack}}},{{_{a}\left( {x,\theta} \right)} = {2{{\ln \left\lbrack \frac{p_{a}\left( {x,\theta} \right)}{p_{a}\left( {\mu_{f},\alpha_{f}} \right)} \right\rbrack}.}}}$

It follows that

l _(ƒ)(x,θ)=M(x,θ;μ,P,α,β,Γ,κ)+(y−h(x,θ))_(T) R⁻¹(y−h(x,θ))−M(μ_(ƒ),α_(ƒ);μ,P,α,β,Γ,κ)−(y−h(μ_(ƒ),α_(ƒ)))^(T) R⁻¹(y−h(μ_(ƒ)α_(ƒ))),

l _(a)(x,θ)=M(x,θ;μ _(ƒ) ,P _(ƒ),α_(ƒ),β_(ƒ),{circumflex over(Γ)}_(ƒ),κ_(ƒ)),

where M is the Mahalanobis von Mises statistic (26). The overdeterminedsystem of algebraic equations

l _(ƒ)(x _(σ) _(i) ,θ_(σ) _(i) )=l _(a)(x _(σ) _(i) ,θ_(σ) _(i) ),i=1, .. . ,N,

is solved for {circumflex over (Γ)}_(ƒ) in the least squares sense. Thisleads to the optimization problem

$\begin{matrix}{\mspace{79mu} {{{\Gamma_{f} = {{\underset{{\hat{\Gamma}}_{f}}{argmin}{\sum\limits_{i = 1}^{n}{\left\lbrack {r_{i}\left( {\hat{\Gamma}}_{f} \right)} \right\rbrack^{2}\mspace{14mu} {subject}\mspace{14mu} {to}\mspace{14mu} {\hat{\Gamma}}_{f}^{T}}}} = {\hat{\Gamma}}_{f}}},\mspace{20mu} {where}}{{{r_{i}\left( {\hat{\Gamma}}_{f} \right)} = {{\xi_{\sigma_{i}}^{T}\xi_{\sigma_{i}}} + {4\kappa \; \sin^{2}\frac{1}{2}\eta_{\sigma_{i}}} + {v_{\sigma_{i}}^{T}v_{\sigma_{i}}} - {\xi_{f}^{T}\xi_{f}} - {4\kappa \; \sin^{2}\frac{1}{2}\eta_{f}} - {v_{f}^{T}v_{f}} - {z_{\sigma_{i}}^{T}z_{\sigma_{i}}} - {4\kappa_{f}\sin^{2}\frac{1}{2}\varphi_{\sigma_{i}}}}},\mspace{20mu} {{{and}\xi_{\sigma_{i}}} = {S^{- 1}\left( {{h\left( {x_{\sigma_{i}},\theta_{\sigma_{i}}} \right)} - y} \right)}},{\xi_{f} = {S^{- 1}\left( {{h\left( {\mu_{f},\alpha_{f}} \right)} - y} \right)}},{R = {SS}^{T}},{v_{\sigma_{i}} = {{\theta_{\sigma_{i}} - {\Theta \left( x_{\sigma_{i}} \right)}} = {\varphi_{\sigma_{i}} + \sigma_{f} + {\beta_{f}^{T}z_{\sigma_{i}}} + {\frac{1}{2}z_{\sigma_{i}}^{T}} - \left( {\alpha + {\beta^{T}v_{\sigma_{i}}} + {\frac{1}{2}v_{\sigma_{i}}^{T}\Gamma \; v_{\sigma_{i}}}} \right)}}},\mspace{20mu} {\eta_{f} = {\alpha_{f} - {\left( {\alpha + {\beta^{T}v_{f}} + {\frac{1}{2}v_{f}^{T}\Gamma \; v_{f}}} \right).}}}}}} & (83)\end{matrix}$

It is noted that only the first two terms in the residuals r_(i), namelyξ_(σ) _(i) ^(T)ξ_(σ) _(i) and 4κ sin² 1/2η_(σ) _(i) , depend on{circumflex over (Γ)}_(ƒ). The other terms can be pre-computed once. Insolving the non-linear least squares problem (83), the startingiteration {circumflex over (Γ)}_(ƒ)=0 is used. It is also useful to notethe partial derivatives of r_(i):

$\frac{\partial r_{i}}{\partial{\hat{\Gamma}}_{f}} = {\left( {{\xi_{\sigma_{i}}^{T}S^{- 1}{\partial_{\theta}{h\left( {x_{\sigma_{i}},\theta_{\sigma_{i}}} \right)}}} + {\kappa \; \sin \; \eta_{\sigma_{i}}}} \right)z_{\sigma_{i}}{z_{\sigma_{i}}^{T}.}}$

At step 304, using the GVM distribution with parameter set (μ_(ƒ),P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ)) computed from Steps 301-303, theprediction error c given by (80) can be expressed in the form

$\begin{matrix}{{{c = {k{\int_{{\mathbb{R}}^{n}}{\int_{- \pi}^{\pi}{\left( {x,{\theta;\mu_{f}},P_{f},\alpha_{f},\beta_{f},\Gamma_{f},\kappa_{f}} \right){f\left( {x,\theta} \right)}{\theta}{x}}}}}},\mspace{20mu} {where}}\mspace{20mu} {k = {\frac{2\pi \; ^{- \kappa_{f}}{I_{0}\left( \kappa_{f} \right)}}{2{\pi }^{- \kappa}{I_{0}(\kappa)}}.\left\lbrack \frac{\det \left( {2\pi \; P_{f}} \right)}{{\det \left( {2\pi \; P} \right)}{\det \left( {2\pi \; R} \right)}} \right\rbrack^{1/2}}}\mspace{20mu} {and}{{{f\left( {x,\theta} \right)} = {\exp \left\lbrack {{{- \frac{1}{2}}z^{T}z} + {\frac{1}{2}z_{f}^{T}z_{f}} - {\frac{1}{2}\xi^{T}\xi} - {2\kappa \; \sin^{2}\frac{1}{2}\left( {\theta - {\Theta (x)}} \right)} + {2\kappa \; \sin^{2}\frac{1}{2}\left( {\theta - {\Theta_{f}(x)}} \right)}} \right\rbrack}},{z = {A^{- 1}\left( {x - \mu} \right)}},{z_{f} = {A_{f}^{- 1}\left( {x - \mu_{f}} \right)}},{\xi = {S^{- 1}\left( {{h\left( {x,\theta} \right)} - y} \right)}},{R = {SS}^{T}},\mspace{20mu} {{\Theta (x)} = {\alpha + {\beta^{T}z} + {\frac{1}{2}z^{T}\Gamma \; z}}},{{\Theta_{f}(x)} = {\alpha_{f} + {\beta_{f}^{T}z_{f}} + {\frac{1}{2}z_{f}^{T}\Gamma_{f}{z_{f}.}}}}}} & (84)\end{matrix}$

The framework of Gauss von Mises quadrature, as developed in Section IV,is applicable to the evaluation of the integral in (84). If {(z_(σ) _(i),φ_(σ) _(i) )}_(i=1) ^(N) and {w_(σ) _(i) }_(i=1) ^(N) are thequadrature nodes and weights for the canonical GVM distribution, then itfollows that

$\begin{matrix}{{{c \approx {k{\sum\limits_{i = 1}^{N}{w_{\sigma_{i}}f_{\sigma_{i}}}}}},{where}}{{f_{\sigma_{i}} = {\exp \begin{bmatrix}{{{- \frac{1}{2}}v_{\sigma_{i}}^{T}v_{\sigma_{i}}} + {\frac{1}{2}z_{\sigma_{i}}^{T}z_{\sigma_{i}}} - {\frac{1}{2}\xi_{\sigma_{i}}^{T}\xi_{\sigma_{i}}} -} \\{{2\kappa \; \sin^{2}\frac{1}{2}\eta_{\sigma_{i}}} + {2\kappa_{f}\sin^{2}\frac{1}{2}\varphi_{\sigma_{i}}}}\end{bmatrix}}},}} & (85)\end{matrix}$

and the definitions of ν_(σ) _(i) , ξ_(σ) _(i) , and η_(σ) _(i) areprovided in Step 303 following Equation (83).

As noted in the case of full state update, the evaluation of theprediction error integral (84) using different order GVM quadraturerules allows the operator to assess when the fused PDF (79) is notwell-represented by a GVM distribution with the parameter set (μ_(ƒ),P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ)). Indeed, the approximation (85) tothe prediction error is exact if the fused PDF (79) is identically a GVMdistribution. In such a case, the function ƒ(x,θ) in (84) can be aconstant. Hence, the prediction error integral can be evaluated exactlyusing any GVM quadrature rule. For example, using a trivial first-orderGVM quadrature rule leads to

${c_{1} = {k\; {\exp \left\lbrack {{{- \frac{1}{2}}v_{f}^{T}v_{f}} - {\frac{1}{2}\xi_{f}^{T}\xi_{f}} - {2\kappa \; \sin^{2}\frac{1}{2}\eta_{f}}} \right\rbrack}}},$

where ν_(ƒ), ξ_(ƒ), and η_(ƒ) are the auxiliary quantities definedfollowing Equation (83). (The subscript ‘1’ underneath c implies the useof a first-order GVM quadrature rule.) Similarly, one can evaluate (84)using a third-order quadrature method in conjunction with Equations(85). The (absolute or relative) difference between the high andlow-order approximations to the prediction error (84) can be used toestimate the error of the integration and validate the underlyingassumptions of the algorithm. A breakdown can be declared if thisdifference exceeds a user-defined threshold.

Stated another way, updating a predicted location of an object in amulti-dimensional space having a plurality of (n+1) dimensions cancomprise reading a first probability density function representing afirst location of the object on a cylindrical manifold (

^(n)×

) of the multi-dimensional space. A second probability density functionrepresenting a second location of the object on the cylindrical manifoldof the multi-dimensional space can be received. The first probabilitydensity function can be fused with the second probability densityfunction to form a third probability density function representing athird location of the object on the cylindrical manifold of themulti-dimensional space. The third probability density function can beprovided as an indication of the updated predicted location of theobject in the multi-dimensional space.

More specifically, each of the first, second, and third probabilitydensity functions can be represented by Gauss von Mises distributionsdefined on the manifold

^(n)×

. The first Gauss von Mises distribution can be represented by theparameter set (μ₁, P₁, α₁, β₁, Γ₁, κ₁), the second Gauss von Misesdistribution can be represented by the parameter set (μ₂, P₂, α₂, β₂,Γ₂, κ₂), and the third Gauss von Mises distribution can be representedby the parameter set (μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ)). Theparameters μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ) can be calculated201 from analytical algebraic and trigonometric expressions. In somecases, a prediction error can be computed 202 using a Gauss von Misesquadrature rule of a chosen order of accuracy. In a SSA application, atleast one of the first or the second Gauss von Mises distributions canbe generated from pluralities of radar, electro-optical, or infraredsensor observations. In various cases, at least one of the first or thesecond Gauss von Mises distributions can be generated from Gaussiandistributions defined on manifold

^(n+1). The second probability density function can comprise anobservation. In such cases, the observation can be related to the firstprobability density function by a stochastic measurement model, and eachof the probability density functions can be defined on a n+1 dimensionalcylindrical manifold (

^(n)×

).

In one embodiment, the first and third probability density functions canbe represented by Gauss von Mises distributions defined on the manifold

^(n)×

. In such cases, the first probability density function can berepresented by the parameter set (μ, P, α, β, Γ, κ) and the thirdprobability density function can be represented by the parameter set(μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ)). The second probabilitydensity function can be a vector (yε

^(p)) hypothesized to be a realization of a random vector (h(x,θ)+v),where h:

^(n)×

→

^(p) and v is a zero-mean p-dimensional Gaussian random vector withcovariance matrix R. The parameters μ_(ƒ) and α_(ƒ) can be computed 301by solving a non-linear least squares problem and the parameters P_(ƒ),β_(ƒ), and κ_(ƒ) can be computed 302 from analytical algebraicexpressions. The parameter Γ_(ƒ) can be computed 303 by solving anon-linear least squares problem and, in some cases, a prediction errorcan be computed 304 using a Gauss von Mises quadrature rule of a chosenorder of accuracy. In a SSA application, the first Gauss von Misesdistribution can be generated from a plurality of radar,electro-optical, or infrared sensor observations. In some cases, thefirst Gauss von Mises distribution can be generated from a Gaussiandistribution defined on the manifold

^(n+1).

VII. Batch Processing

Sequential non-linear filtering, the focus of the previous two sections,entails updating the system state recursively using multiple uncertaintypropagation (prediction) and data fusion (correction) operations withina Bayesian framework. This section considers the batch filtering problemwhich processes a sequence of reports simultaneously to produce aprobability density function (PDF) of the system state conditioned onthe reports. Such batch processing capabilities can be used for trackinitiation or orbit determination to produce a prior PDF of the systemstate which can be subsequently passed to a sequential filter.

For a system state defined on a cylindrical manifold, such as a spaceobject's equinoctial orbital element state, a Bayesian maximum aposteriori framework can be used to generate a Gauss von Mises (GVM)distribution representing the state conditioned on a sequence ofreports. The input reports can be either (i) GVM distributions in thesame state space or (ii) observations each related to the system stateby a stochastic measurement model. When specialized to two-body dynamicsin orbital mechanics, the new batch processing capability using the GVMdistribution can extend the sequential differential correction (SEQDC)and batch differential correction (BATCHDC) algorithms in theAstrodynamics Standards software suite maintained by Air Force SpaceCommand. SEQDC and BATCHDC are the analogous algorithms for Cases (i)and (ii), respectively, and produce a Gaussian PDF of the orbital stateconditioned on the reports. The complete algorithm descriptions forCases (i) and (ii) are provided in Subsections VII.A and VII.B,respectively, with respective block diagrams provided in FIGS. 4 and 5.

VII.A. Full State Reports

FIG. 4 is a flowchart illustrating an exemplary process for generating aGauss von Mises distribution from a plurality of reports, wherein saidreports are Gauss von Mises distributions according to one embodiment ofthe present invention. In the case of full state reports, the batchprocessing algorithm can have the following inputs 400:

-   -   i. a sequence of mutually independent random vectors (x_(k),        θ_(k)) ε        ^(n)×        , for k=1, . . . , m, called the reports, each distributed as a        GVM distribution with respective parameter set (μ_(k), P_(k),        α_(k), β_(k), Γ_(k), κ_(k)); and    -   ii. a sequence of diffeomorphisms Φ_(k):        ^(n)×        →        ^(n)×        , for k=1, . . . , m−1.        The output 404 can be a random vector (x_(m), θ_(m)) ε        ^(n)×        , distributed as a GVM distribution with parameter set ({tilde        over (μ)}_(m), {tilde over (P)}_(m), {tilde over (α)}_(m),        {tilde over (β)}_(m), {tilde over (Γ)}_(m), {tilde over        (κ)}_(m)) which characterizes the “fusion” of the m-th report        with the first m−1 reports subject to the constraints

(x _(k+1),θ_(k+1))=Φ_(k)(x _(k),θ_(k)),

for k=1, . . . , m−1. The algorithm can have three main steps summarizedbelow:

-   -   Compute 401 the components {tilde over (μ)}_(m) and {tilde over        (α)}_(m) of the output GVM parameter set by solving the least        squares problem (90);    -   Recover 402 the components {tilde over (P)}_(m), {tilde over        (β)}_(m), and {tilde over (κ)}_(m) of the output GVM parameter        set from Equation (91); and    -   Compute 403 the component {tilde over (Γ)}_(m) of the output GVM        parameter set by solving the least squares problem (93). If it        is believed that the output is well-approximated by a Gaussian,        this step can be skipped and one can set {tilde over (Γ)}_(m)=0.        These three steps are described in more detail below including        discussions on theoretical considerations and a precise        statistical definition of what is meant by “fusion” in the sense        of the output.

VII.A.i. Theoretical Considerations

In a general setting, suppose for k=1, . . . , m that X_(k) is a randomvector defined on a manifold

with PDF p_(X) _(k) _(|I) _(k) (x_(k)|I_(k)), where I_(k) denotes anyprior information association with X_(k). Let

_(m)=(I₁, . . . , I_(m)). The joint PDF of X₁, . . . , X_(m)(conditioned on

_(m)) can be denoted by

p _(X) ₁ _(, . . . ,X) _(m) _(|)

_(m) (x ₁ , . . . ,x _(m)|

_(m)).

The conditional PDF of X_(m) given X_(k)=x_(k), for k=1, . . . , m−1,and given the prior information

_(m) can be denoted by

p _(X) _(m) _(|X) ₁ _(, . . . ,X) _(m−1) _(,)

_(m) (x _(m) |x ₁ , . . . ,x _(m−1),

_(m)).

Let Φ_(k):

→

, for k=1, . . . , m−1, be a sequence of diffeomorphisms withcorresponding actions

x _(k+1)=Φ_(k)(x _(k)),k=1, . . . ,m−1,

and let Ψ_(k) denote the inverse of Φ_(k) so that

x _(k)=Ψ_(k)(x _(k+1)),k=1, . . . ,m−1.

Further, for k=1, . . . , m−1, define the composition maps

Ξ_(m−k)=Ψ_(m−k)∘Ψ_(m−k+1)∘ . . . ∘Ψ_(m−1)  (86)

with corresponding actions

x _(m−k)=Ξ_(m−k)(x _(m)),k=1, . . . ,m−1

For notational convenience, Ξ_(m) denotes the identity map.

The PDF obtained by fusing the m-th report with the first m−1 reportscan be defined to be the conditional PDF of X_(m) givenX_(k)=Ξ_(k)(x_(m)), for k=1, . . . , m−1, and given the other priorinformation

_(m):

p _(ƒ)(x _(m|)

_(m))=p _(X) _(m) _(|X) ₁ _(, . . . ,X) _(m−1)

_(m)(x _(m)|Ξ₁(x _(m)), . . . ,Ξ_(m−1)(x _(m)),

_(m)).

By the definition of conditional probability,

p f  ( x m | m ) = 1 c  p X 1 ,  …  , X m | m  ( Ξ 1  ( x m ) , … , Ξ m  ( x m ) | m ) , ( 87 )

where the normalization constant c is given by

c =  p X 1 ,  …  , X m - 1 | m  ( Ξ 1  ( x m ) , …  , Ξ m - 1  (x m ) | m ) =  ∫ ℳ  p X 1 ,  …  , X m | m  ( Ξ 1  ( x m ) , …  ,Ξ m  ( x m ) | m )    x m . ( 88 )

Additionally, if (i) the sequence of random vectors X₁, . . . , X_(m)are mutually independent given the prior information

_(m) and (ii) for k=1, . . . , m, X_(k) is independent of I₁, . . . ,I_(k−1), I_(k+1), . . . , I_(m) given I_(k), then (87) and (88) simplifyto

p f  ( x m | m ) = 1 c  ∏ k = 1 m   p X k | I k  ( Ξ k  ( x m ) |I k ) ,  c = ∫ ℳ   ∏ k = 1 m  p X k | I k  ( Ξ k  ( x m ) | I k )  x m .

The above theoretical considerations and independence assumptions cannow be specialized to case when the random vectors X₁, . . . , X_(m) arethe input described at the beginning of the subsection (and I₁, . . . ,I_(m) denote any prior information used to construct their respectiveparameter sets).

VII.A.ii. Algorithm Description

At step 401, given the input defined at the beginning of the subsection,the fused PDF in (x_(m),θ_(m)) is

$\begin{matrix}{{{p_{f}\left( {x_{m},\theta_{m}} \right)} = {\frac{1}{c}{\prod\limits_{k = 1}^{m}\; {{\mathcal{M}}\left( {{{\Xi_{k}\left( {x_{m},\theta_{m}} \right)};\mu_{k}},P_{k},\alpha_{k},\beta_{k},\Gamma_{k},\kappa_{k}} \right)}}}},} & (89)\end{matrix}$

where c is a normalization constant and the Ξ_(k), for k=1, . . . , m,are the composition maps defined in (86). The objective can be toapproximate (89) by a GVM distribution with parameters ({tilde over(μ)}_(m), {tilde over (P)}_(m), {tilde over (α)}_(m), {tilde over(β)}_(m), {tilde over (Γ)}_(m), {tilde over (κ)}_(m)). The mode of thisoutput GVM distribution is ({tilde over (μ)}_(m), {tilde over (α)}_(m)).These parameters can be chosen to coincide with the mode of (89). Thecomputation of the mode can be formulated in terms of a non-linear leastsquares problem as follows. Define

${{{q_{f}\left( {x_{m},\theta_{m}} \right)} \equiv {{- \ln}\; {p_{f}\left( {x_{m},\theta_{m}} \right)}}} = {{\frac{1}{2}{r\left( {x_{m},\theta_{m}} \right)}^{T}{r\left( {x_{m},\theta_{m}} \right)}} + k}},$

where k is a constant (independent of x_(m) and θ_(m)), and

${{r\left( {x_{m},\theta_{m}} \right)} = \begin{bmatrix}z_{1} \\{2\sqrt{\kappa_{1}}\sin \; \frac{1}{2}\varphi_{1}} \\\vdots \\z_{m} \\{2\sqrt{\kappa_{m}}\sin \; \frac{1}{2}\varphi_{m}}\end{bmatrix}},{where}$${z_{k} = {A_{k}^{- 1}\left( {x_{k} - \mu_{k}} \right)}},{\varphi_{k} = {\theta_{k} - \alpha_{k} - {\beta_{k}^{T}z_{k}} - {\frac{1}{2}z_{k}^{T}\Gamma_{k}z_{k}}}},{\left( {x_{k},\theta_{k}} \right) = {\Xi_{k}\left( {x_{m},\theta_{m}} \right)}},$

for i=1, . . . , m. Thus, the parameters μ_(ƒ) and α_(ƒ) are computed bysolving the non-linear least squares problem

$\begin{matrix}\begin{matrix}{\left( {{\overset{\sim}{\mu}}_{m},{\overset{\sim}{\alpha}}_{m}} \right) = {\begin{matrix}{argmax} \\{x_{m},\theta_{m}}\end{matrix}{p_{f}\left( {x_{m},\theta_{m}} \right)}}} \\{= {\begin{matrix}{argmin} \\{x_{m},\theta_{m}}\end{matrix}\frac{1}{2}{r\left( {x_{m},\theta_{m}} \right)}^{T}{{r\left( {x_{m},\theta_{m}} \right)}.}}}\end{matrix} & (90)\end{matrix}$

At step 402, the parameters {tilde over (P)}_(m), {tilde over (β)}_(m),and {tilde over (κ)}_(m) of the output GVM distribution can be recoveredfrom the osculating Gaussian of (89) in analogy to Step 302 of the datafusion algorithm in Section VI.B. Specifically,

$\begin{matrix}{{\begin{bmatrix}{{{\overset{\sim}{A}}_{m}^{- T}\left( {I + {{\overset{\sim}{\kappa}}_{m}{\overset{\sim}{\beta}}_{m}{\overset{\sim}{\beta}}_{m}^{T}}} \right)}A_{m}^{- 1}} & {{- {\overset{\sim}{\kappa}}_{m}}{\overset{\sim}{A}}_{m}^{- T}{\overset{\sim}{\beta}}_{m}} \\{{- \kappa_{m}}{\overset{\sim}{\beta}}_{m}^{T}{\overset{\sim}{A}}_{m}^{- 1}} & {\overset{\sim}{\kappa}}_{m}\end{bmatrix} = {\partial_{x_{m}}^{2}{q_{f}\left( {{\overset{\sim}{\mu}}_{m},{\overset{\sim}{\alpha}}_{m}} \right)}}},} & (91)\end{matrix}$

where {tilde over (P)}_(m)=Ã_(m)Ã_(m) ^(T) and

_(m)=(x_(m), θ_(m)). Analytic expressions for the second derivatives ofq_(ƒ) in terms of the derivatives of the maps Ξ_(k) are given asfollows. First, the notation Ξ_(kx), Ξ_(kθ), and Ξ_(ki), for i=1, . . ., n, is used to express the map Ξ_(k) in component form according to

${{\Xi_{k}\left( {x_{m},\theta_{m}} \right)} = {\begin{bmatrix}{\Xi_{kx}\left( {x_{m},\theta_{m}} \right)} \\{\Xi_{k\; \theta}\left( {x_{m},\theta_{m}} \right.}\end{bmatrix} = \begin{bmatrix}{\Xi_{k\; 1}\left( {x_{m},\theta_{m}} \right)} \\\vdots \\{\Xi_{kn}\left( {x_{m},\theta_{m}} \right)} \\{\Xi_{k\; \theta}\left( {x_{m},\theta_{m}} \right)}\end{bmatrix}}},$

for k=1, . . . , m. Next, define the auxiliary quantities

B _(k) =A _(k) ⁻¹∂

_(m) Ξ_(kx) ,b _(k)=∂

_(m) Ξ_(kθ) −B _(k) ^(T)(β_(k)+Γ_(k) z _(k)),

c _(k) =A _(k) ^(−T) z _(k) ,d _(k) =A _(k) ^(−T)(β_(k)+Γ_(k) z _(k)),

for k=1, . . . , m. The second derivatives of q_(ƒ) are given by

${{\partial_{x_{m}}^{2}q_{f}} = {{\left( {\partial_{x_{m}}r} \right)^{T}{\partial_{x_{m}}r}} + {\sum\limits_{i = 1}^{{({n + 1})}m}{r_{i}{\partial_{x_{m}}^{2}r_{i}}}}}},$

where r_(i) denotes the i-th component of the residual r. The Jacobianof r with respect to

_(m) is

${\partial_{x_{m}}r} = {\begin{bmatrix}B_{1} \\b_{1}^{T} \\\vdots \\B_{m} \\b_{M}^{T}\end{bmatrix}.}$

The second derivative terms can be decomposed according to

${\sum\limits_{i = 1}^{{({n + 1})}m}{r_{i}{\partial_{x_{m}}^{2}r_{i}}}} = {{\sum\limits_{k = 1}^{m}{\sum\limits_{i = 1}^{n}{\left( z_{k} \right)_{i}{\partial_{x_{m}}^{2}\left( z_{k} \right)_{i}}}}} + {\sum\limits_{k = 1}^{m}{2\sqrt{\kappa_{k}}\sin \; \frac{1}{2}\varphi_{k}{{\partial_{x_{m}}^{2}\left( {2\sqrt{\kappa_{k}}\sin \; \frac{1}{2}\varphi_{k}} \right)}.}}}}$

It follows that

$\mspace{20mu} {{{\sum\limits_{i = 1}^{n}{\left( z_{k} \right)_{i}{\partial_{x_{m}}^{2}\left( z_{k} \right)_{i}}}} = {\sum\limits_{i = 1}^{n}{\left( c_{k} \right)_{i}{\partial_{x_{m}}^{2}\Xi_{k\; i}}}}},\mspace{20mu} {and}}$${\partial_{x_{m}}^{2}\left( {2\sqrt{\kappa_{k}}\sin \; \frac{1}{2}\varphi_{k}} \right)} = {{\sqrt{\kappa_{k}}\left\lbrack {{\cos \; \frac{1}{2}{\varphi_{k}\left( {{\partial_{x_{m}}^{2}\Xi_{k\; \theta}} - {B_{k}^{T}\Gamma_{k}B_{k}} - {\sum\limits_{i = 1}^{n}{\left( d_{k} \right)_{i}{\partial_{x_{m}}^{2}\Xi_{k\; i}}}}} \right)}} - {\frac{1}{2}\sin \; \frac{1}{2}\varphi_{k}b_{k}b_{k}^{T}}} \right\rbrack}.}$

These partial derivatives are also useful if solving the optimizationproblem (90) in Step 401 using Gauss-Newton or full Newton iterations.

At step 403, let Γ _(m) be an approximation to {tilde over (Γ)}_(m) andinitially set T _(m)=0. In this step, the approximation Γ _(m) can berefined by solving a non-linear least squares problem, in analogy toStep 104 of the uncertainty propagation algorithm in Section V.C andStep 303 of the data fusion algorithm in Section VI.B. As motivated inthe latter, one can elect to skip this step (and set {tilde over(Γ)}_(m)=0) if it is believed that the output is well-approximated by aGaussian, as might be the case if the underlying application is thespace surveillance tracking problem.

Using the exact expression for the fused PDF p_(ƒ)(x_(m), θ_(m)) givenby (89) and the approximation given by

p _(a)(x _(m),θ_(m))=

(x _(m),θ_(m);{tilde over (μ)}_(m) ,{tilde over (P)} _(m),{tilde over(α)}_(m),{tilde over (β)}_(m),{tilde over (Γ)}_(m),{tilde over(κ)}_(m)),  (92)

the approximation of Γ _(m) is improved by “matching” p_(a) and p_(ƒ) atthe N quadrature nodes {(x_(σ) _(i) ,θ_(σ) _(i) )}_(i=1) ^(N) of the GVMdistribution with parameter set ({tilde over (μ)}_(m), {tilde over(P)}_(m), {tilde over (α)}_(m), {tilde over (β)}_(m), Γ _(m), {tildeover (κ)}_(m)). These quadrature nodes are defined from either thethird, fifth, or higher-order GVM quadrature rules derived in SectionsIV.A, IV.B, and IV.C, respectively. It is noted that

${x_{\sigma_{i}} = {{\overset{\sim}{\mu}}_{m} + {{\overset{\sim}{A}}_{m}z_{\sigma_{i}}}}},{\theta_{\sigma_{i}} = {\varphi_{\sigma_{i}} + {\overset{\sim}{\alpha}}_{m} + {{\overset{\sim}{\beta}}_{m}^{T}z_{\sigma_{i}}} + {\frac{1}{2}z_{\sigma_{i}}^{T}{\overset{\_}{\Gamma}}_{m}z_{\sigma_{i}}}}},$

for i=1, . . . , N, where {(z_(σ) _(i) ,φ_(σ) _(i) )}_(i=1) ^(N) are thequadrature nodes corresponding to the canonical GVM distribution. Oneway to perform this matching procedure is to consider the matchingequations in ‘log space’ by defining

${{_{f}\left( {x_{m},\theta_{m}} \right)} = {{- 2}\; {\ln \left\lbrack \frac{p_{f}\left( {x_{m},\theta_{m}} \right)}{p_{f}\left( {{\overset{\sim}{\mu}}_{m},{\overset{\sim}{\alpha}}_{m}} \right)} \right\rbrack}}},{{_{a}\left( {x_{m},\theta_{m}} \right)} = {{- 2}{{\ln \left\lbrack \frac{p_{a}\left( {x_{m},\theta_{m}} \right)}{p_{a}\left( {{\overset{\sim}{\mu}}_{m},{\overset{\sim}{\alpha}}_{m}} \right)} \right\rbrack}.}}}$

It follows that

${{_{f}\left( {x_{m},\theta_{m}} \right)} = {\sum\limits_{k = 1}^{m}\begin{bmatrix}{{M\left( {{{\Xi_{k}\left( {x_{m},\theta_{m}} \right)};\mu_{k}},P_{k},\alpha_{k},\beta_{k},\Gamma_{k},\kappa_{k}} \right)} -} \\{M\left( {{{\Xi_{k}\left( {{\overset{\sim}{\mu}}_{m},{\overset{\sim}{\alpha}}_{m}} \right)};\mu_{k}},P_{k},\alpha_{k},\beta_{k},\Gamma_{k},\kappa_{k}} \right)}\end{bmatrix}}},\mspace{20mu} {{_{a}\left( {x_{m},\theta_{m}} \right)} = {M\left( {x_{m},{\theta_{m};{\overset{\sim}{\mu}}_{m}},{\overset{\sim}{P}}_{m},{\overset{\sim}{\alpha}}_{m},{\overset{\sim}{\beta}}_{m},{\overset{\sim}{\Gamma}}_{m},{\overset{\sim}{\kappa}}_{m}} \right)}},$

where M is the Mahalanobis von Mises statistic (26). The overdeterminedsystem of algebraic equations

l _(ƒ)(x _(σ) _(i) ,θ_(σ) _(i) )=l _(a)(x _(σ) _(i) ,θ_(σ) _(i) ),i=1, .. . ,N,

is solved for Γ _(m) in the least squares sense. This leads to theoptimization problem

$\begin{matrix}{\mspace{79mu} {{{{\overset{\sim}{\Gamma}}_{m} = {{\underset{{\overset{\_}{\Gamma}}_{m}}{argmin}{\sum\limits_{i = 1}^{n}{\left\lbrack {r_{i}\left( {\overset{\_}{\Gamma}}_{m} \right)} \right\rbrack^{2}\mspace{14mu} {subject}\mspace{14mu} {to}\mspace{14mu} {\overset{\_}{\Gamma}}_{m}^{T}}}} = {\overset{\_}{\Gamma}}_{m}}},\mspace{79mu} {where}}{{{r_{i}\left( {\overset{\_}{\Gamma}}_{m} \right)} = {{\sum\limits_{k = 1}^{m}\left\lbrack {{v_{k,\sigma_{i}}^{T}v_{k,\sigma_{i}}} + {4\; \kappa_{k}\sin^{2}\frac{1}{2}\eta_{k,\sigma_{i}}} - {{\overset{\sim}{v}}_{k}^{T}{\overset{\sim}{v}}_{k}} - {4\; \kappa_{k}\sin^{2}\frac{1}{2}{\overset{\sim}{\eta}}_{k}}} \right\rbrack} - {z_{\sigma_{i}}^{T}z_{\sigma_{i}}} - {4\; {\overset{\sim}{\kappa}}_{m}\sin^{2}\frac{1}{2}\varphi_{\sigma_{i}}}}},\mspace{79mu} {and}}{{v_{k,\sigma_{i}} = {A_{k}^{- 1}\left( {{\Xi_{k\mspace{14mu} x}\left( {x_{\sigma_{i}},\theta_{\sigma_{i}}} \right)} - \mu_{k}} \right)}},{{\overset{\sim}{v}}_{k} = {A_{k}^{- 1}\left( {{\Xi_{k\mspace{14mu} x}\left( {{\overset{\sim}{\mu}}_{m},{\overset{\sim}{\alpha}}_{m}} \right)} - \mu_{k}} \right)}},\mspace{20mu} {\eta_{k,\sigma_{i}} = {{\Xi_{k\mspace{14mu} \theta}\left( {x_{\sigma_{i}},\theta_{\sigma_{i}}} \right)} - \alpha_{k} - {\beta_{k}^{T}v_{k,\sigma_{i}}} - {\frac{1}{2}v_{k,\sigma_{i}}^{T}\Gamma_{k}v_{k,\sigma_{i}}}}},\mspace{20mu} {{\overset{\sim}{\eta}}_{k} = {{\Xi_{k\mspace{14mu} \theta}\left( {{\overset{\sim}{\mu}}_{m},{\overset{\sim}{\alpha}}_{m}} \right)} - \alpha_{k} - {\beta_{k}^{T}{\overset{\sim}{v}}_{k}} - {\frac{1}{2}{\overset{\sim}{v}}_{k}^{T}\Gamma_{k}{{\overset{\sim}{v}}_{k}.}}}}}}} & (93)\end{matrix}$

It is noted that the only terms in the residuals r_(i) that depend on Γ_(m) are ν_(k,σ) _(i) ^(T)ν_(k,σ) _(i) and

$4\; \kappa_{k}\sin^{2}\frac{1}{2}{\eta_{k,\sigma_{i}}.}$

The other terms can be pre-computed once. In solving the non-linearleast squares problem (93), the starting iteration Γ _(m)=0 is used. Itis also useful to note the partial derivatives of r_(i):

$\frac{\partial r_{i}}{\partial{\overset{\_}{\Gamma}}_{m}} = {z_{\sigma_{i}}z_{\sigma_{i}}^{T}{\sum\limits_{k = 1}^{m}{\begin{Bmatrix}{{v_{k,\sigma_{i}}A_{k}^{- 1}{\partial_{\theta}{\Xi_{k\mspace{14mu} x}\left( {x_{\sigma_{i}},\theta_{\sigma_{i}}} \right)}}} + \quad} \\{\kappa_{k}\sin \; \eta_{k,\sigma_{i}} \times \begin{bmatrix}{{\partial_{\theta}{\Xi_{k\mspace{14mu} \theta}\left( {x_{\sigma_{i}},\theta_{\sigma_{i}}} \right)}} -} \\{\left( {\beta_{k} + {\Gamma_{k}z_{\sigma_{i}}}} \right)^{T}A_{k}^{- 1}{\partial_{\theta}{\Xi_{k\mspace{14mu} x}\left( {x_{\sigma_{i}},\theta_{\sigma_{i}}} \right)}}}\end{bmatrix}}\end{Bmatrix}.}}}$

Optionally, the fidelity of the approximation (92) can be verified byevaluating the normalization constant c using a first and third-orderaccurate GVM quadrature rule, in analogy to the discussions at the endof the data fusion algorithms in Sections VIA and VI.B. A breakdown isdeclared if the (absolute or relative) difference exceeds a user-definedthreshold.

VII.B. Measurement-Based Reports

FIG. 5 is a flowchart illustrating an exemplary process for generating aGauss von Mises distribution from a plurality of reports, wherein saidreports are observations related to the state space by a stochasticmeasurement model according to one embodiment of the present invention.In the case of measurement-based reports, the batch processing algorithmcan have the following inputs 500:

-   -   i. a sequence of smooth maps h_(k):        ^(n)×        →        ^(p), called the measurement models,    -   ii. a sequence of random vectors v_(k)ε        ^(p), for k=1, . . . , m, called the measurement errors,        distributed jointly as a zero-mean Gaussian random vector of        dimension pm with covariance        ,    -   iii. a sequence of diffeomorphisms Φ_(k):        ^(n)×        →        ^(n)×        , for k=1, . . . , m−1,    -   iv. a sequence of vectors y_(k)ε        ^(p), for k=1, . . . , m, called the measurements or reports,        each hypothesized to be a realization of the respective random        vector h_(k)(x_(k),θ_(k))+v_(k), where the (x_(k),θ_(k)) are        constrained by the relations        (x_(k+1),θ_(k+1))=Φ_(k)(x_(k),θ_(k)), for k=1, . . . , m−1, and    -   v. optionally, a random vector ({circumflex over (x)}_(m),        {circumflex over (θ)}_(m)) ε        ^(n)×        , called the prior, distributed as a GVM distribution with        parameter set ({circumflex over (μ)}_(m), {circumflex over        (P)}_(m), {circumflex over (α)}_(m), {circumflex over (β)}_(m),        {circumflex over (Γ)}_(m), {circumflex over (κ)}_(m)).        The output 502 can be a random vector (x_(m),θ_(m)) ε        ^(n)×        , distributed as a GVM distribution with parameter set (μ_(m),        P_(m), α_(m), β_(m), Γ_(m), κ_(m)) which characterizes the        “fusion” of the m reports and the optional prior.

In what follows, theoretical considerations are discussed and a precisestatistical definition of what is meant by “fusion” in the sense of theoutput is made. In the process, it is shown that batch processing withmeasurement-based reports in conjunction with the input listed abovereduces to the data fusion algorithm of Section VI.B for ameasurement-based update. This is Step 501 in FIG. 5.

VII.B.i. Theoretical Considerations

In a general setting, suppose X_(m) is a random vector, called theprior, defined on a manifold

with PDF p_(X) _(m) _(|I) _(m) (x_(m)|I_(m)), where I_(m) denotes anyprior information associated with X_(m). Further, suppose V₁, . . . ,V_(k) are random vectors each defined on a manifold

with joint PDF p_(V) ₁ _(, . . . ,V) _(m) (v₁, . . . , v_(m)). Given asequence of smooth maps h_(k):

→

, for k=1, . . . , m, and a sequence of diffeomorphisms Φ_(k):

→

, for k=1, . . . , m−1, define the sequence of random vectors Y_(k), fork=1, . . . , m, according to

Y _(k) =h _(k)(X _(k))+V _(k),

subject to the constraints

X _(k+1)=Φ_(k)(X _(k)),

for k=1, . . . , m−1. Let Ψ_(k) denote the inverse of Φ_(k) so that

X _(k)=Ψ_(k)(X _(k+1)),k=1, . . . ,m−1.

Further, for k=1, . . . , m−1, define the composition maps

Ξ_(m−k)=Ψ_(m−k)∘Ψ_(m−k+1)∘ . . . ∘Ψ_(m−1),  (94)

and let Ξ_(m) be the identity map. For notational convenience, let

=(V₁, . . . , V_(m)),

=(X₁, . . . , X_(m)),

=(x₁, . . . , x_(m)),

=(Y₁, . . . , Y_(m)), and

=(y₁, . . . , y_(m)).

The PDF obtained by fusing a sequence of reports y₁, . . . , y_(m)ε

with the prior is defined to be the conditional PDF of X_(m) given thatY_(k)=y_(k), for k=1, . . . , m, and given the other prior informationI_(m):

p _(ƒ)(x _(m)|

,I_(m))≡p _(X) _(m) _(|Y,I) _(m) (x _(m)|

,I_(m)).

By marginalization,

p _(ƒ)(x _(m)|

,I_(m))=

. . .

p_(X|Y,I) _(m) (

|

,I_(m))dx ₁ . . . dx _(m−1),  (95)

By Bayes' rule,

${{p_{{x|},I_{m}}\left( {\left. x \right|,I_{m}} \right)} = {\frac{1}{c}{p_{{|x},I_{m}}\left( {\left| x \right.,I_{m}} \right)}{p_{x|I_{m}}\left( x \middle| I_{m} \right)}}},$

where c=p_(Y|I) _(m) (

|I_(m)). From the preceding definitions, it follows that

${{p_{{|x},I_{m}}\left( {\left| x \right.,I_{m}} \right)} = {p_{v}\left( {{y_{1} - {h_{1}\left( x_{1} \right)}},\ldots \mspace{14mu},{y_{m} - {h_{m}\left( x_{m} \right)}}} \right)}},{\quad\begin{matrix}{{p_{x|I_{m}}\left( x \middle| I_{m} \right)} = {{p_{X_{m}|I_{m}}\left( x_{m} \middle| I_{m} \right)}{\prod\limits_{k = 1}^{m - 1}{p_{{X_{k}|X_{k + 1}},I_{m}}\left( {\left. x_{k} \middle| x_{k + 1} \right.,I_{m}} \right)}}}} \\{{= {{p_{X_{m}|I_{m}}\left( x_{m} \middle| I_{m} \right)}{\prod\limits_{k = 1}^{m - 1}{\delta \left( {x_{k} - {\Psi_{k}\left( x_{k + 1} \right)}} \right)}}}},}\end{matrix}}$

where ‘δ’ denotes the Dirac delta function. Substituting theseexpressions into (95) yields

${{p_{f}\left( {\left. x_{m} \right|,I_{m}} \right)} = {\frac{1}{c}{p_{v}\left( {{y_{1} - {\left( {h_{1}\bullet \mspace{11mu} \Xi_{1}} \right)\left( x_{m} \right)}},\ldots \mspace{14mu},{y_{m} - {\left( {h_{m}\bullet \mspace{11mu} \Xi_{m}} \right)\left( x_{m} \right)}}} \right)}{p_{X_{m}|I_{m}}\left( x_{m} \middle| I_{m} \right)}}},$

where the normalization constant c is given by

c=

p _(V)(y ₁−(h ₁∘Ξ₁)(x _(m)), . . . ,y _(m)−(h _(m)∘Ξ_(m))(x _(m)))p _(X)_(m) _(|I) _(m) (x _(m) |I _(m))dx _(m).

If there is no prior information on X_(m), then the prior PDF p_(X) _(m)_(|I) _(m) (x_(m)|I_(m)) is removed from the formulation by making thediffuse prior assumption.

The above theoretical considerations are now specialized to the casewhen the random vectors X₁, . . . , X_(m) and V₁, . . . , V_(m), themaps h₁, . . . , h_(m) and Φ₁, . . . , Φ_(m−1), and the reports y₁, . .. , y_(m) are the input described at the beginning of the subsection. Itis shown how the algorithm of Section VI.B can be applied to generatethe desired output.

VII.B.ii. Algorithm Description

At step 501, given the input defined at the beginning of the subsection,the fused PDF in (x_(m),θ_(m)) is

$\begin{matrix}{{{p_{f}\left( {x_{m},\theta_{m}} \right)} = {\frac{1}{c}\left( {{{- {\left( {x_{m},\theta_{m}} \right)}};0},} \right)\left( {x_{m},{\theta_{m};{\hat{\mu}}_{m}},{\hat{P}}_{m},{\hat{\alpha}}_{m},{\hat{\beta}}_{m},{\hat{\Gamma}}_{m},{\hat{\kappa}}_{m}} \right)}},} & (96)\end{matrix}$

where c is a normalization constant and

$\begin{matrix}{{= \begin{bmatrix}y_{1} \\\vdots \\y_{m}\end{bmatrix}},} & {= {\begin{bmatrix}{h_{1}\bullet \mspace{11mu} \Xi_{1}} \\\vdots \\{h_{m}\bullet \mspace{11mu} \Xi_{m}}\end{bmatrix}.}}\end{matrix}$

The composition maps Ξ_(k), for k=1, . . . , m, can be defined in termsof the input diffeomorphisms Φ₁, . . . , Φ_(m−1) by way of Equations(86).

The objective is to approximate (96) by a GVM distribution withparameter set (μ_(m), P_(m), α_(m), β_(m), Γ_(m), κ_(m)). This isprecisely the objective of the data fusion algorithm in Section VI.B fora measurement-based update upon comparing Equations (96) and (79).Indeed, the algorithm of Section VI.B is applicable in this case wherethe inputs are (i) the prior GVM distribution with parameter set({circumflex over (μ)}_(m), {circumflex over (P)}_(m), {circumflex over(α)}_(m), {circumflex over (β)}_(m), {circumflex over (Γ)}_(m),{circumflex over (κ)}_(m)), (ii) the map

defined above, (iii) the covariance matrix

, and (iv) the vector

defined above. If there is no prior information on the state (x_(m),θ_(m)) for the first input, then any terms dependent on the prior in theobjective function of Equation (81) can be removed.

Stated another way, determining a track of an object through amulti-dimensional space having a plurality of (n+1) dimensions cancomprise receiving a plurality of probability density functionsincluding at least one report probability density function. Theplurality of probability density functions can be fused to form anoutput probability density function. Each probability density functioncan be defined on a n+1 dimensional cylindrical manifold (

^(n)×

), and the probability density functions can be constrained by aplurality of diffeomophisms from

^(n)×

to

^(n)×

. The output probability density function can be provided as arepresentation of the track of the object through the multi-dimensionalspace.

For example, each of the plurality of probability density functions cancomprise a report represented by a Gauss von Mises distribution definedon the manifold

^(n)×

and the Gauss von Mises distributions can comprise at least a firstdistribution with a parameter sets (μ_(k), P_(k), α_(k), β_(k), Γ_(k),κ_(k)), and a second distribution with a parameter set ({tilde over(μ)}_(m), {tilde over (P)}_(m), {tilde over (α)}_(m), {tilde over(β)}_(m), {tilde over (Γ)}_(m), {tilde over (κ)}_(m)). In such cases,the parameters {tilde over (μ)}_(m) and {tilde over (α)}_(m) be computed401 by solving a non-linear least squares problem. The parameters {tildeover (P)}_(m), {tilde over (β)}_(m), and {tilde over (κ)}_(m) can becomputed 402 from analytical algebraic expressions. The parameter {tildeover (Γ)}_(m) can be computed 403 by solving a non-linear least squaresproblem. In a SSA application, one or more of the Gauss von Misesdistributions can be generated from pluralities of radar,electro-optical, or infrared sensor observations, one or more of theGauss von Mises distributions can be generated from Gaussiandistributions defined on the manifold

^(n+1). The diffeomorphisms can be solution flows induced from a systemof ordinary differential equations describing two-body dynamics oforbital mechanics for the object.

In another case, the plurality of probability density functions cancomprise a first probability density function and at least one report.Each probability density function can be represented by a Gauss vonMises distribution defined on the manifold

^(n)×

. The first Gauss von Mises distribution can be represented by aparameter set comprising ({circumflex over (μ)}_(m), {circumflex over(P)}_(m), {circumflex over (α)}_(m), {circumflex over (β)}_(m),{circumflex over (Γ)}_(m), {circumflex over (κ)}_(m)) and the reportGauss von Mises distributions can be represented by a parameter setcomprising (μ_(m), P_(m), α_(m), β_(m), Γ_(m), κ_(m)). The report Gaussvon Mises distributions can be vectors (y_(k)ε

^(p), for k=1, . . . , m), each hypothesized to be a realization of arandom vector (h(x_(k),θ_(k))+v_(k)), where h_(k):

^(n)×

→

^(p), for k=1, . . . , m, and (v₁, . . . , v_(m)) is a zero-meanpm-dimensional Gaussian random vector with covariance matrix

. In a SSA application, the first Gauss von Mises distribution can begenerated from a plurality of radar, electro-optical, or infrared sensorobservations. The first Gauss von Mises distribution can be generatedfrom a Gaussian distribution defined on the manifold

^(n+1). One or more of the reports may also be generated from radar,electro-optical, or infrared sensors. The diffeomorphisms can besolution flows induced from a system of ordinary differential equationsdescribing two-body dynamics of orbital mechanics for the object.

EXAMPLES

In the following, an example is presented to demonstrate embodiments ofthe present invention. The example is a scenario in low Earth orbit(LEO) which compares the Gauss von Mises (GVM) filter prediction stepwith that of the extended Kalman filter (EKF), the unscented Kalmanfilter (UKF), a Gaussian sum filter (GSF), and a particle filter. Theaccuracy of the GVM uncertainty propagation algorithm is also validatedusing a metric based on the L₂ error. For the specific LEO scenario, itis shown that the GVM filter prediction step properly characterizes theactual uncertainty of a space object's orbital state while simple lesssophisticated methods which make Gaussian assumptions (such as the EKFand UKF) do not. Specifically, under the non-linear propagation oftwo-body dynamics, the new algorithm properly characterizes theuncertainty for up to eight times as long as the standard UKF all at noadditional computational cost to the UKF.

In what follows, the particulars of the simulation scenario are definedin Section I and the results are discussed Section II.

I. Scenario Description

This section describes the specific input to the uncertainty propagationalgorithm in Section V.C above, how the output is visualized, and howthe accuracy of the output is validated.

I.A. Input

The initial GVM distribution of the space object's orbital state (i.e.,input (i)) is defined with respect to equinoctial orbital elements (a,h, k, p, q, l)ε

⁵×

with parameter set (μ, P, α, β, Γ, κ) given by

$\begin{matrix}{{\mu = \begin{bmatrix}{7136.635\mspace{14mu} {km}} \\0 \\0 \\0 \\0\end{bmatrix}},} & {P = \begin{bmatrix}\left( {20\mspace{14mu} {km}} \right)^{2} & 0 & 0 & 0 & 0 \\0 & 10^{- 6} & 0 & 0 & 0 \\0 & 0 & 10^{- 6} & 0 & 0 \\0 & 0 & 0 & 10^{- 6} & 0 \\0 & 0 & 0 & 0 & 10^{- 6}\end{bmatrix}}\end{matrix},\mspace{20mu} {\alpha = 0},{\beta = 0},{\Gamma = 0},{\kappa = {3.282806 \times {10^{7}.}}}$

The mode of this distribution describes a circular, non-inclined orbitin LEO with a semi-major axis of 7136.635 km. This choice of semi-majoraxis is made so that the instantaneous orbital period of the object is100 minutes. In subsequent discussions, a time unit of one orbitalperiod is equal to 100 minutes. It is noted that the GVM distributionwith the above parameter set is approximately Gaussian (since Γ=0 andκ>>1). In particular, the standard deviation in the mean longitudecoordinate l is σ_(l)=1/√{square root over (κ)}=0.01°=36″. Whenvalidating against the EKF, UKF, or Gaussian sum filter, the osculatingGaussian (22) is used to convert the input GVM distribution to aGaussian distribution.

The diffeomorphism describing the non-linear transformation (i.e., input(ii)) is the solution flow (47) induced from the system of ordinarydifferential equations (2) describing the two-body dynamics of orbitalmechanics. (See the discussions in Sections I.A and V.A above) Theparameters of the diffeomorphism are the epoch time t₀ of the input, thefinal epoch time t, and the specific forces used to model theperturbations. In simulations, t₀ is the J2000 epoch (1 Jan. 2000, 12:00UTC) and t−t₀ is varied from 0.5 to 8 orbital periods. The EGM96 gravitymodel of degree and order 70 is used to model the perturbations. Thoughone can include additional non-conservative forces such as atmosphericdrag, the non-linearities induced by the gravity alone (especially inLEO) are sufficiently strong to stress the algorithm. Moreover, theinitial standard deviation in the semi-major axis of 20 km ispessimistic. It is observed that non-Gaussian effects are acceleratedfor larger initial uncertainties in the semi-major axis. Finally, thenumerical integration of the ordinary differential equations isperformed using a Gauss-Jackson method.

I.B. Visualization

The output of the uncertainty propagation algorithm in Section V.C.above is a GVM distribution characterizing the uncertainty of the spaceobject's orbital state at a specified future epoch. In order tovisualize this six-dimensional probability density function (PDF), thelevel curves (i.e., curves of equal likelihood) of the two-dimensional(2D) marginal PDF in the semi-major axis a and mean longitudecoordinates l are plotted. This choice is made because it is along thisparticular 2D slice where the greatest departure from “Gaussianity” andthe most extreme “banana” or “boomerang” shaped level curves areobserved.

In the panels of FIG. 9, nσ level curves (of the marginal PDF in the aand l coordinates) are plotted in half-sigma increments starting at

$n = \frac{1}{2}$

for various final epoch times t. In order to visualize these marginalPDFs which are defined on a cylinder, the cylinder is cut at {tilde over(α)}−π (where {tilde over (α)} is the α parameter of the propagated GVMdistribution) and rolled out to form a 2D plane. This plane is thenrotated so that the semi-major and semi-minor axes of the osculatingGaussian covariance are aligned with the horizontal and vertical. (Anysuch rotation or rescaling does not exaggerate any non-Gaussian effectsor the extremity of the boomerang shape because cumulants of order threeand higher are invariant under an affine transformation.) Whereappropriate, overlays of the level curves generated from the EKF and UKFare also shown. Additionally, the crosses represent particles generatedfrom a Monte-Carlo based particle filter. If the represented PDFproperly characterizes the actual uncertainty, then approximately 99.5%of the particles should be contained within the respective 3σ levelcurve.

I.C. Validation

An inspection of the panels in FIG. 9 provides a simple visual means toassess whether the represented uncertainty properly characterizes theactual uncertainty of the space object's orbital state. “Most” of theparticles (crosses) should be contained within the level curves. Toquantify uncertainty realism more rigorously and hence validate theprediction steps of the different filters under consideration, thenormalized L₂ error is studied.

For functions ƒ,g:

→

, the normalized L₂ error between ƒ and g is the scalar

${{{\overset{\_}{L}}_{2}\left( {f,g} \right)} = \frac{{{f - g}}_{L_{2}}^{2}}{{f}_{L_{2}}^{2} + {g}_{L_{2}}^{2}}},$

where ∥•∥_(L) ₂ is the L₂ norm:

∥ƒ∥_(L) ₂ ²=

ƒ(x)² dx.

By non-negativity of the L₂ norm, L ₂≧0 with equality if ƒ=g in the L₂sense. By the triangle inequality, L ₂≦1 with equality if ƒ and g areorthogonal in the L₂ sense.

The validation tests shown in FIG. 10 generate a time history of thenormalized L₂ error

L ₂(p _(approx)(•,t),p _(baseline)(•,t)),

where t is the final epoch time. Further, p_(approx)(u,t) represents anapproximation to the PDF of a space object's orbital state at time t (urepresents equinoctial orbital elements) as computed by the predictionstep of the GVM filter, EKF, UKF, or a Gaussian sum filter. Moreover,p_(baseline)(u,t) is a high-fidelity approximation to the exact statePDF which serves as a baseline for comparison. This baseline is computedusing a high-fidelity GVM mixture distribution using the Gaussian sumrefinement scheme of Horwood et al. extended to GVM distributions asdescribed in Section V.F above. Thus, papProx is a good approximation tothe actual state PDF if the normalized L₂ error between it and thebaseline Pbaseiine is zero or “close to zero” for all time.

II. Discussions

As described in Subsection I.B. of the EXAMPLE, the panels in FIGS. 9a-9 f show the evolution of a space object's orbital uncertainty (withinitial conditions defined in Subsection I.A of the EXAMPLE) computedusing the prediction steps of the EKF, UKF, GVM filter, and a particlefilter. Each of the six panels of FIGS. 9 a-9 f shows the respectivelevel curves 925 for EKF, 920 for UKF, and 915 for GVM in the plane ofthe semi-major axis and mean longitude coordinates at the epochs t−t₀=0,0.5, 1, 2, 4, 8 orbital periods. In each of the six epochs, the levelcurves 915 produced by the GVM filter correctly capture the actualuncertainty depicted by the particle ensemble. Each set of level curves915, 920, and 925 is deduced from the propagation of only 13 quadraturenodes (corresponding to a third-order GVM quadrature rule). Thecomputational cost is the same as that of the UKF. For the UKF, itscovariance is indeed consistent (realistic) in the sense that it agreeswith that computed from the definition of the covariance. Thus, in thisscenario, the UKF provides “covariance realism” but clearly does notsupport “uncertainty realism” since the covariance does not representthe actual banana-shaped uncertainty of the exact PDF. Further, thestate estimate produced from the UKF coincides with the mean of theexact PDF. However, the mean is displaced from the mode. Consequently,the probability that the object is within a small neighborhood centeredat the UKF state estimate (mean) is essentially zero. The EKF, on theother hand, provides a state estimate coinciding closely with the mode,but the covariance tends to collapse making inflation necessary to beginto cover the uncertainty. In neither the EKF nor UKF case does thecovariance actually model the uncertainty. In summary, the GVM filtermaintains a proper characterization of the uncertainty, the EKF and UKFdo not.

FIG. 10 shows the evolution of the normalized L₂ error, as defined inSubsection I.C of the EXAMPLE, computed from the UKF 1025, EKF 1020, theGaussian sum filter (GSF) with N=17 and N=49 components 1015 and 1005,and the GVM filter 1010. The UKF and EKF quickly break down, butaccuracy can be improved by increasing the fidelity of the Gaussian sum.The normalized L₂ errors produced from the GVM filter lie between thoseproduced from the 17 and 49-term Gaussian sum. It is noted that the17-term Gaussian sum requires the propagation of 17×13=221 quadraturenodes. The GVM filter uses 13. In principle, the accuracy of the GVMfilter could be further improved by using mixtures of GVM distributionsas described in Section V.F above. Notwithstanding these comments, ifone deems a normalized L₂ error of L ₂=0.05 to signal a breakdown inaccuracy, then the UKF and EKF first hit this threshold after about oneorbital period. By examining when the normalized L₂ error first crosses0.05 for the GVM filter, it is seen that one can propagate theuncertainty using the GVM filter for about eight times longer than whenusing an EKF or UKF.

FIG. 11 is a block diagram illustrating components of an exemplaryoperating environment in which various embodiments of the presentinvention may be implemented. The system 1100 can include one or moreuser computers 1105, 1110, which may be used to operate a client,whether a dedicate application, web browser, etc. The user computers1105, 1110 can be general purpose personal computers (including, merelyby way of example, personal computers and/or laptop computers runningvarious versions of Microsoft Corp.'s Windows and/or Apple Corp.'sMacintosh operating systems) and/or workstation computers running any ofa variety of commercially-available UNIX or UNIX-like operating systems(including without limitation, the variety of GNU/Linux operatingsystems). These user computers 1105, 1110 may also have any of a varietyof applications, including one or more development systems, databaseclient and/or server applications, and web browser applications.Alternatively, the user computers 1105, 1110 may be any other electronicdevice, such as a thin-client computer, Internet-enabled mobiletelephone, and/or personal digital assistant, capable of communicatingvia a network (e.g., the network 1115 described below) and/or displayingand navigating web pages or other types of electronic documents.Although the exemplary system 1100 is shown with two user computers, anynumber of user computers may be supported.

In some embodiments, the system 1100 may also include a network 1115.The network can be any type of network familiar to those skilled in theart that can support data communications using any of a variety ofcommercially-available protocols, including without limitation TCP/IP,SNA, IPX, AppleTalk, and the like. Merely by way of example, the network1115 may be a local area network (“LAN”), such as an Ethernet network, aToken-Ring network and/or the like; a wide-area network; a virtualnetwork, including without limitation a virtual private network (“VPN”);the Internet; an intranet; an extranet; a public switched telephonenetwork (“PSTN”); an infra-red network; a wireless network (e.g., anetwork operating under any of the IEEE 802.11 suite of protocols, theBluetooth protocol known in the art, and/or any other wirelessprotocol); and/or any combination of these and/or other networks such asGSM, GPRS, EDGE, UMTS, 3G, 2.5 G, CDMA, CDMA2000, WCDMA, EVDO etc.

The system may also include one or more server computers 1120, 1125,1130 which can be general purpose computers and/or specialized servercomputers (including, merely by way of example, PC servers, UNIXservers, mid-range servers, mainframe computers rack-mounted servers,etc.). One or more of the servers (e.g., 1130) may be dedicated torunning applications, such as a business application, a web server,application server, etc. Such servers may be used to process requestsfrom user computers 1105, 1110. The applications can also include anynumber of applications for controlling access to resources of theservers 1120, 1125, 1130.

The web server can be running an operating system including any of thosediscussed above, as well as any commercially-available server operatingsystems. The web server can also run any of a variety of serverapplications and/or mid-tier applications, including HTTP servers, FTPservers, CGI servers, database servers, Java servers, businessapplications, and the like. The server(s) also may be one or morecomputers which can be capable of executing programs or scripts inresponse to the user computers 1105, 1110. As one example, a server mayexecute one or more web applications. The web application may beimplemented as one or more scripts or programs written in anyprogramming language, such as Java™, C, C# or C++, and/or any scriptinglanguage, such as Perl, Python, or TCL, as well as combinations of anyprogramming/scripting languages. The server(s) may also include databaseservers, including without limitation those commercially available fromOracle®, Microsoft®, Sybase®, IBM® and the like, which can processrequests from database clients running on a user computer 1105, 1110.

In some embodiments, an application server may create web pagesdynamically for displaying on an end-user (client) system. The web pagescreated by the web application server may be forwarded to a usercomputer 1105 via a web server. Similarly, the web server can receiveweb page requests and/or input data from a user computer and can forwardthe web page requests and/or input data to an application and/or adatabase server. Those skilled in the art will recognize that thefunctions described with respect to various types of servers may beperformed by a single server and/or a plurality of specialized servers,depending on implementation-specific needs and parameters.

The system 1100 may also include one or more databases 1135. Thedatabase(s) 1135 may reside in a variety of locations. By way ofexample, a database 1135 may reside on a storage medium local to (and/orresident in) one or more of the computers 1105, 1110, 1115, 1125, 1130.Alternatively, it may be remote from any or all of the computers 1105,1110, 1115, 1125, 1130, and/or in communication (e.g., via the network1120) with one or more of these. In a particular set of embodiments, thedatabase 1135 may reside in a storage-area network (“SAN”) familiar tothose skilled in the art. Similarly, any necessary files for performingthe functions attributed to the computers 1105, 1110, 1115, 1125, 1130may be stored locally on the respective computer and/or remotely, asappropriate. In one set of embodiments, the database 1135 may be arelational database, such as Oracle 10 g, that is adapted to store,update, and retrieve data in response to SQL-formatted commands.

FIG. 12 illustrates an exemplary computer system 1200, in which variousembodiments of the present invention may be implemented. The system 1200may be used to implement any of the computer systems described above.The computer system 1200 is shown comprising hardware elements that maybe electrically coupled via a bus 1255. The hardware elements mayinclude one or more central processing units (CPUs) 1205, one or moreinput devices 1210 (e.g., a mouse, a keyboard, etc.), and one or moreoutput devices 1215 (e.g., a display device, a printer, etc.). Thecomputer system 1200 may also include one or more storage device 1220.By way of example, storage device(s) 1220 may be disk drives, opticalstorage devices, solid-state storage device such as a random accessmemory (“RAM”) and/or a read-only memory (“ROM”), which can beprogrammable, flash-updateable and/or the like.

The computer system 1200 may additionally include a computer-readablestorage media reader 1225 a, a communications system 1230 (e.g., amodem, a network card (wireless or wired), an infra-red communicationdevice, etc.), and working memory 1240, which may include RAM and ROMdevices as described above. In some embodiments, the computer system1200 may also include a processing acceleration unit 1235, which caninclude a DSP, a special-purpose processor and/or the like.

The computer-readable storage media reader 1225 a can further beconnected to a computer-readable storage medium 1225 b, together (and,optionally, in combination with storage device(s) 1220) comprehensivelyrepresenting remote, local, fixed, and/or removable storage devices plusstorage media for temporarily and/or more permanently containingcomputer-readable information. The communications system 1230 may permitdata to be exchanged with the network 1220 and/or any other computerdescribed above with respect to the system 1200.

The computer system 1200 may also comprise software elements, shown asbeing currently located within a working memory 1240, including anoperating system 1245 and/or other code 1250, such as an applicationprogram (which may be a client application, web browser, mid-tierapplication, RDBMS, etc.). It should be appreciated that alternateembodiments of a computer system 1200 may have numerous variations fromthat described above. For example, customized hardware might also beused and/or particular elements might be implemented in hardware,software (including portable software, such as applets), or both.Further, connection to other computing devices such as networkinput/output devices may be employed. Software of computer system 1200may include code 1250 for implementing embodiments of the presentinvention as described herein.

In the foregoing description, for the purposes of illustration, methodswere described in a particular order. It should be appreciated that inalternate embodiments, the methods may be performed in a different orderthan that described. It should also be appreciated that the methodsdescribed above may be performed by hardware components or may beembodied in sequences of machine-executable instructions, which may beused to cause a machine, such as a general-purpose or special-purposeprocessor or logic circuits programmed with the instructions to performthe methods. These machine-executable instructions may be stored on oneor more machine readable mediums, such as CD-ROMs or other type ofoptical disks, floppy diskettes, ROMs, RAMs, EPROMs, EEPROMs, magneticor optical cards, flash memory, or other types of machine-readablemediums suitable for storing electronic instructions. Alternatively, themethods may be performed by a combination of hardware and software.

While illustrative and presently preferred embodiments of the inventionhave been described in detail herein, it is to be understood that theinventive concepts may be otherwise variously embodied and employed, andthat the appended claims are intended to be construed to include suchvariations, except as limited by the prior art.

What is claimed is:
 1. A method for updating a predicted location of anobject in a multi-dimensional space having a plurality of (n+1)dimensions, the method comprising: reading, by a computer system, afirst probability density function representing a first location of theobject on a cylindrical multi-dimensional space (

^(n)×

) of the multi-dimensional space; receiving, by the computer system, asecond probability density function representing a second location ofthe object on the cylindrical multi-dimensional space of themulti-dimensional space; fusing, by the computer system, the firstprobability density function with the second probability densityfunction to form a third probability density function representing athird location of the object on the cylindrical multi-dimensional spaceof the multi-dimensional space; and providing, by the computer system,the third probability density function as an indication of the updatedpredicted location of the object in the multi-dimensional space.
 2. Themethod of claim 1, wherein each of the first, second, and thirdprobability density functions are represented by Gauss von Misesdistributions defined on the multi-dimensional space

^(n)×

, wherein the first Gauss von Mises distribution is represented by theparameter set (μ₁, P₁, α₁, β₁, Γ₁, κ₁), wherein the second Gauss vonMises distribution is represented by the parameter set (μ₂, P₂, α₂, β₂,Γ₂, κ₂), and wherein the third Gauss von Mises distribution isrepresented by the parameter set (μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ),κ_(ƒ)).
 3. The method of claim 2, further comprising computing, by thecomputer system, the parameters μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ)from analytical algebraic and trigonometric expressions.
 4. The methodof claim 2, further comprising computing, by the computer system, aprediction error using a Gauss von Mises quadrature rule of a chosenorder of accuracy.
 5. The method of claim 2, wherein at least one of thefirst or the second Gauss von Mises distributions are generated frompluralities of radar, electro-optical, or infrared sensor observations.6. The method of claim 2, wherein at least one of the first or thesecond Gauss von Mises distributions are generated from Gaussiandistributions defined on multi-dimensional space

^(n+1).
 7. The method of claim 1, wherein the second probability densityfunction comprises an observation, wherein the observation is related tothe first probability density function by a stochastic measurementmodel, and wherein each of the probability density functions are definedon a n+1 dimensional cylindrical multi-dimensional space (

^(n)×

).
 8. The method of claim 7, wherein the first and third probabilitydensity functions are represented by Gauss von Mises distributionsdefined on the multi-dimensional space

^(n)×

, wherein the first probability density function is represented by theparameter set (μ, P, α, β, Γ, κ) and the third probability densityfunction is represented by the parameter set (μ_(ƒ), P_(ƒ), α_(ƒ),β_(ƒ), Γ_(ƒ), κ_(ƒ)), and wherein said second probability densityfunction is a vector (yε

^(p)) hypothesized to be a realization of a random vector (h(x,θ)+v),wherein h:

^(n)×

→

^(p) and v is a zero-mean p-dimensional Gaussian random vector withcovariance matrix R.
 9. The method of claim 8, further comprisingcomputing, by the computer system, the parameters μ_(ƒ) and α_(ƒ) bysolving a non-linear least squares problem.
 10. The method of claim 8,further comprising computing, by the computer system, the parametersP_(ƒ), β_(ƒ), and κ_(ƒ) from analytical algebraic expressions.
 11. Themethod of claim 8, further comprising computing, by the computer system,the parameter Γ_(ƒ) by solving a non-linear least squares problem. 12.The method of claim 8, further comprising computing, by the computersystem, the prediction error using a Gauss von Mises quadrature rule ofa chosen order of accuracy.
 13. The method of claim 8, wherein the firstGauss von Mises distribution is generated from a plurality of radar,electro-optical, or infrared sensor observations.
 14. The method ofclaim 8, wherein the first Gauss von Mises distribution is generatedfrom a Gaussian distribution defined on the multi-dimensional space

^(n+1).
 15. The method of claim 8, wherein the observation is generatedfrom a radar, electro-optical, or infrared sensor.
 16. A systemcomprising: a processor; and a memory coupled with and readable by theprocessor and having stored therein a sequence of instruction which,when executed by the processor, cause the processor to update apredicted location of an object in a multi-dimensional space having aplurality of (n+1) dimensions by: reading a first probability densityfunction representing a first location of the object on a cylindricalmulti-dimensional space (

^(n)×

) of the multi-dimensional space, receiving a second probability densityfunction representing a second location of the object on the cylindricalmulti-dimensional space of the multi-dimensional space, fusing the firstprobability density function with the second probability densityfunction to form a third probability density function representing athird location of the object on the cylindrical multi-dimensional spaceof the multi-dimensional space, and providing the third probabilitydensity function as an indication of the updated predicted location ofthe object in the multi-dimensional space.
 17. The system of claim 16,wherein each of the first, second, and third probability densityfunctions are represented by Gauss von Mises distributions defined onthe multi-dimensional space

^(n)×

, wherein the first Gauss von Mises distribution is represented by theparameter set (μ₁, P₁, α₁, β₁, Γ₁, κ₁), wherein the second Gauss vonMises distribution is represented by the parameter set (μ₂, P₂, α₂, β₂,Γ₂, κ₂), and wherein the third Gauss von Mises distribution isrepresented by the parameter set (μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ),κ_(ƒ)).
 18. The system of claim 17, further comprising computing theparameters μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ) from analyticalalgebraic and trigonometric expressions.
 19. The system of claim 17,further comprising computing a prediction error using a Gauss von Misesquadrature rule of a chosen order of accuracy.
 20. The system of claim17, wherein at least one of the first or the second Gauss von Misesdistributions are generated from pluralities of radar, electro-optical,or infrared sensor observations.
 21. The system of claim 17, wherein atleast one of the first or the second Gauss von Mises distributions aregenerated from Gaussian distributions defined on multi-dimensional space

^(n+1).
 22. The system of claim 16, wherein the second probabilitydensity function comprises an observation, wherein the observation isrelated to the first probability density function by a stochasticmeasurement model, and wherein each of the probability density functionsare defined on a n+1 dimensional cylindrical multi-dimensional space (

^(n)×

).
 23. The system of claim 22, wherein the first and third probabilitydensity functions are represented by Gauss von Mises distributionsdefined on the multi-dimensional space

^(n)×

, wherein the first probability density function is represented by theparameter set (μ, P, α, β, Γ, κ) and the third probability densityfunction is represented by the parameter set (μ_(ƒ), P_(ƒ), α_(ƒ),β_(ƒ), Γ_(ƒ), κ_(ƒ)), and wherein said second probability densityfunction is a vector (yε

^(p)) hypothesized to be a realization of a random vector (h(x,θ)+v),wherein h:

^(n)×

→

^(p) and v is a zero-mean p-dimensional Gaussian random vector withcovariance matrix R.
 24. The system of claim 23, further comprisingcomputing the parameters μ_(ƒ) and α_(ƒ) by solving a non-linear leastsquares problem.
 25. The system of claim 23, further comprisingcomputing the parameters P_(ƒ), β_(ƒ), and κ_(ƒ) from analyticalalgebraic expressions.
 26. The system of claim 23, further comprisingcomputing the parameter Γ_(ƒ)by solving a non-linear least squaresproblem.
 27. The system of claim 23, further comprising computing theprediction error using a Gauss von Mises quadrature rule of a chosenorder of accuracy.
 28. The system of claim 23, wherein the first Gaussvon Mises distribution is generated from a plurality of radar,electro-optical, or infrared sensor observations.
 29. The system ofclaim 23, wherein the first Gauss von Mises distribution is generatedfrom a Gaussian distribution defined on the multi-dimensional space

^(n+1).
 30. The system of claim 23, wherein the observation is generatedfrom a radar, electro-optical, or infrared sensor.
 31. Acomputer-readable memory having stored therein a sequence ofinstructions which, when executed by a processor, cause the processor toupdate a predicted location of an object in a multi-dimensional spacehaving a plurality of (n+1) dimensions by: reading a first probabilitydensity function representing a first location of the object on acylindrical multi-dimensional space (

^(n)×

) of the multi-dimensional space; receiving a second probability densityfunction representing a second location of the object on the cylindricalmulti-dimensional space of the multi-dimensional space; fusing the firstprobability density function with the second probability densityfunction to form a third probability density function representing athird location of the object on the cylindrical multi-dimensional spaceof the multi-dimensional space; providing the third probability densityfunction as an indication of the updated predicted location of theobject in the multi-dimensional space.
 32. The computer-readable memoryof claim 31, wherein each of the first, second, and third probabilitydensity functions are represented by Gauss von Mises distributionsdefined on the multi-dimensional space

^(n)×

, wherein the first Gauss von Mises distribution is represented by theparameter set (μ₁, P₁, α₁, β₁, Γ₁, κ₁), wherein the second Gauss vonMises distribution is represented by the parameter set (μ₂, P₂, α₂, β₂,Γ₂, κ₂), and wherein the third Gauss von Mises distribution isrepresented by the parameter set (μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ),κ_(ƒ)).
 33. The computer-readable memory of claim 32, further comprisingcomputing the parameters μ_(ƒ), P_(ƒ), α_(ƒ), β_(ƒ), Γ_(ƒ), κ_(ƒ) fromanalytical algebraic and trigonometric expressions.
 34. Thecomputer-readable memory of claim 32, further comprising computing aprediction error using a Gauss von Mises quadrature rule of a chosenorder of accuracy.
 35. The computer-readable memory of claim 32, whereinat least one of the first or the second Gauss von Mises distributionsare generated from pluralities of radar, electro-optical, or infraredsensor observations.
 36. The computer-readable memory of claim 32,wherein at least one of the first or the second Gauss von Misesdistributions are generated from Gaussian distributions defined onmulti-dimensional space

^(n+)1.
 37. The computer-readable memory of claim 31, wherein the secondprobability density function comprises an observation, wherein theobservation is related to the first probability density function by astochastic measurement model, and wherein each of the probabilitydensity functions are defined on a n+1 dimensional cylindricalmulti-dimensional space (

^(n)×

).
 38. The computer-readable memory of claim 37, wherein the first andthird probability density functions are represented by Gauss von Misesdistributions defined on the multi-dimensional space

^(n)×

, wherein the first probability density function is represented by theparameter set (μ, P, α, β, Γ, κ) and the third probability densityfunction is represented by the parameter set (μ_(ƒ), P_(ƒ), α_(ƒ),β_(ƒ), Γ_(ƒ), κ_(ƒ)), and wherein said second probability densityfunction is a vector (yε

) hypothesized to be a realization of a random vector (h(x,θ)+v),wherein h:

^(n)×

→

and v is a zero-mean p-dimensional Gaussian random vector withcovariance matrix R.
 39. The computer-readable memory of claim 38,further comprising computing the parameters μ_(ƒ) and α_(ƒ) by solving anon-linear least squares problem.
 40. The computer-readable memory ofclaim 38, further comprising computing the parameters P_(ƒ), β_(ƒ), andκ_(ƒ) from analytical algebraic expressions.
 41. The computer-readablememory of claim 38, further comprising computing the parameter Γ_(ƒ) bysolving a non-linear least squares problem.
 42. The computer-readablememory of claim 38, further comprising computing the prediction errorusing a Gauss von Mises quadrature rule of a chosen order of accuracy.43. The computer-readable memory of claim 38, wherein the first Gaussvon Mises distribution is generated from a plurality of radar,electro-optical, or infrared sensor observations.
 44. Thecomputer-readable memory of claim 38, wherein the first Gauss von Misesdistribution is generated from a Gaussian distribution defined on themulti-dimensional space

^(n+1).
 45. The computer-readable memory of claim 38, wherein theobservation is generated from a radar, electro-optical, or infraredsensor.